Experimental validation of computerised models of clustering of platelet glycoprotein receptors that signal via tandem SH2 domain proteins

The clustering of platelet glycoprotein receptors with cytosolic YxxL and YxxM motifs, including GPVI, CLEC-2 and PEAR1, triggers activation via phosphorylation of the conserved tyrosine residues and recruitment of the tandem SH2 (Src homology 2) domain effector proteins, Syk and PI 3-kinase. We have modelled the clustering of these receptors with monovalent, divalent and tetravalent soluble ligands and with transmembrane ligands based on the law of mass action using ordinary differential equations and agent-based modelling. The models were experimentally evaluated in platelets and transfected cell lines using monovalent and multivalent ligands, including novel nanobody-based divalent and tetravalent ligands, by fluorescence correlation spectroscopy. Ligand valency, receptor number, receptor dimerisation, receptor phosphorylation and a cytosolic tandem SH2 domain protein act in synergy to drive receptor clustering. Threshold concentrations of a CLEC-2-blocking antibody and Syk inhibitor act in synergy to block platelet aggregation. This offers a strategy for countering the effect of avidity of multivalent ligands and in limiting off-target effects.


Introduction
Platelets are regulated by two major classes of surface receptors: seven transmembrane receptors, which signal through heterotrimeric G proteins, and glycoprotein receptors, which signal through tyrosine kinases. Current antiplatelet drugs such as aspirin and P2Y 12 receptor antagonists target G protein-coupled receptor (GPCR) regulated pathways but are not effective in all individuals and cause a significant increase in risk of bleeding in some people and this can be life-threatening. The glycoprotein receptors GPVI, CLEC-2 and PEAR1 are considered to be targets for a new class of antithrombotic agents based on mouse models and genome-wide association studies with inhibition predicted to cause less bleeding than current antiplatelet agents [1,2]. These receptors are activated by clustering, and antagonists must therefore overcome the net effect of ligand affinity and avidity for effective inhibition, thus favouring inhibitors with high affinity and a slow off-rate of dissociation. A deeper understanding of the relationship between ligand engagement and receptor clustering of this class of receptor will guide development of more potent antithrombotic agents.
The major signalling glycoprotein receptors in platelets have tyrosine-based signalling motifs in their cytosolic tails which bind to the tandem SH2 domain proteins Syk and phosphatidylinositol (PI) 3-kinase. This includes the GPVI-FcRγ complex and the low affinity immune receptor, FcγRIIA, which signal via an immuno-receptor-tyrosine-based activation motif (ITAM) characterised by two conserved YxxL sequences separated by 6-12 amino acids [3,4]. CLEC-2, on the other hand, signals via a single YxxL known as a hemITAM, and PEAR1 via a single YxxM sequence [2,5]. The binding of the tandem SH2 domains in Syk to phosphorylated ITAM or two hemITAM sequences initiates a signalling cascade that leads to activation of phospholipase C (PLC) γ2. The binding of the tandem SH2 domains in PI 3-kinase to two phosphorylated YxxM sequences leads to the activation of PI 3-kinase and generation of phosphatidylinositol 3,4,5-trisphosphate (PIP3) and recruitment of pleckstrin homology (PH) domain proteins to the membrane, including the kinases Akt and Btk.
In platelets, a net increase in phosphorylation of the cytosolic tails of glycoprotein receptors occurs when a sufficient density is reached for Src and Syk kinases to overcome the effect of active tyrosine phosphatases in the membrane such as CD148, SHP1 and SHP2 [6,7]. Multivalent ligands therefore induce activation through receptor clustering whereas monovalent ligands act as antagonists. In addition, ligand engagement may lead to a conformational change that increases the availability of tyrosines for phosphorylation, as exemplified by the dissociation of ITAMs from the inner leaflet of the membrane in the T cell receptor [8]. There is no evidence however that this plays a role in mediating activation of platelet glycoprotein receptors.
The reversible binding of a monovalent ligand to a receptor has been modelled by the law of mass action which states that the rate of a chemical reaction is proportional to the concentration of the reactants [9]. As originally proposed, the law of mass action describes receptor occupancy but not the conformational change that gives rise to the activation of multi-membrane spanning receptors such as GPCRs and ion channels. This was followed by the two-state model of receptor theory which accounts for this conformational change and provides a molecular explanation of efficacy [10]. More complex models of receptor activation have since emerged to account for additional features of receptors such as desensitisation and inverse agonists (for review see [11]).
In 2008, Cooper and Qian [12] applied the law of mass action to monomeric receptors which signal through Src family kinases. The ordinary differential equations (ODEs) developed in their study led them to conclude that phosphorylation of these receptors by Src kinases can occur solely as a result of receptor clustering without the need for a net change in kinase activity, receptor conformation or movement to a specialised location in the membrane such as a lipid raft. They proposed that clustering of receptors generates a critical density of receptor tyrosines for phosphorylation and binding of Src kinases via their SH2 domains, promoting transphosphorylation of the Src kinases and receptor phosphorylation thereby initiating downstream signalling [12]. However, they also emphasised that additional factors could play a role, such as other signalling proteins or the cytoskeleton. These additional interactions can be further developed using ODE-based models, but they rapidly become complex.
Agent-based modelling (ABM) is an alternative computational method that does not require a full understanding of all of the components of the system and can be used to simulate the actions and interactions of autonomous agents in order to identify key regulatory steps in a complex system [13,14]. ABM also has the advantage that it enables spatial relationships to be considered, whereas this is not easily done for equation-based methods. ABM has been used in many areas, including in sociology, town-planning and economics, but has only been applied to the study of receptors in a few instances, with a primary focus on the modelling of complex intracellular signalling pathways [15,16]. In ABM, it is relatively easy to incorporate rules that govern the clustering of receptors such as ligand binding, location, phosphorylation and interactions with other receptors and signalling molecules, and thereby study their interactions.
In the present study, we have used ODE and ABM to explore the reversible binding of soluble ligands (monovalent, divalent and tetravalent) and membrane-bound ligands (monovalent) to platelet glycoprotein receptors. For the ODE-based models, we have assumed that the receptors are monomers, but for the agent-based models, we have modelled receptors as a mixture of monomers and dimers, and introduced further complexity through a divalent cytosolic crosslinker to mimic the tandem SH2 domains of Syk and PI 3-kinase. We have tested the models by monitoring receptor activation in platelets and by measurement of receptor clustering in transfected cell lines using fluorescence correlation spectroscopy (FCS) in the presence of monovalent, divalent and multivalent ligands.  Fig  1A and 1B. Examples of monovalent ligands are an antibody Fab or a nanobody, and of a divalent ligand an antibody or dimerised nanobody. A tetravalent ligand can be generated by the crosslinking of monovalent ligands such as a nanobody, as generated in this study, and is representative of a multivalent ligand. Fig 1C illustrates how a tandem SH2 domain-containing protein, representing Syk or PI 3-kinase, can further support receptor clustering by causing cross-linkage of receptors in the cytoplasm. Fig 1D and 1E show a monovalent membrane ligand and an immobilised multivalent ligand, as exemplified by the CLEC-2 ligand podoplanin and the GPVI ligand collagen, respectively.

Monomeric receptors can be clustered by a variety of ligands
The law of mass action was introduced to model the interaction of a drug with its receptor by Clark in 1933 [9] with the following assumptions: (i) all reactions are reversible (ii) the interactions are one-to-one (iii) each interaction is an independent event (iv) the ligand concentration is in excess over that of the receptor such that the change in ligand concentration on binding to the receptor is negligible. These assumptions enabled the relationship between ligand concentration and response to be investigated and later, when radioligand binding studies were developed, the link to occupancy. The second and third assumptions are not valid, however, for glycoprotein receptors that are activated by clustering as receptor binding and activation are dependent on the combination of affinity and avidity, while the fourth assumption does not apply to surface-restricted ligands. Nevertheless, the principle that the interaction of a single epitope in a ligand with a monomeric receptor is determined by the law of mass action and that this is the basis of activation has been shown to be valid [12]. The present study further explores this relationship using ODEs and ABM as applied to ITAM, hemI-TAM and YxxM-containing receptors.

ODE modelling
ODE modelling of the interaction of a monovalent ligand and a monomeric receptor. The interaction of a monovalent soluble ligand and a monomeric receptor was modelled using ODEs according to law of mass action [17]. At equilibrium, the law defines the amount of the ligand-receptor [LR] complex in relation to the concentration of ligand [L] and receptor [R] as shown in Equation (Eq) 1.1: The amount of LR is determined by the association (k 1 ) and dissociation (k -1 ) rate constants, which have units of M -1 s -1 and s -1 , respectively. In the model by Clark [9], the concentration of the ligand is in vast excess of the receptor and it follows that, at equilibrium, the system reaches a steady state defined by the equilibrium dissociation constant, K D (Eq 1.2). The unit of K D is M, and K D is given by: This may also be expressed as: where R tot denotes the total number of receptors. The reciprocal of K D is the equilibrium association constant, K A , with units of M -1 : For the model where the ligand is assumed to be in vast excess of the receptor, and therefore the concentration is not significantly reduced, the K D is equivalent to the concentration of ligand required to occupy 50% of the receptors, with greater than 95% receptor occupancy occurring at a 19-fold higher concentration. The dynamics of the ligand-receptor interaction are defined by a linear ODE, whose solution gives the following standard result [18]: The relationship between ligand concentration and occupancy when plotted on a semi-log plot is sigmoid (Fig 2A).
The assumption that the ligand concentration is in excess over the receptor concentration does not apply to a transmembrane ligand since the concentration will reduce upon binding to the receptor, unless it can be rapidly replenished from other regions of the membrane or by fusion of intracellular membranes. Thus, the membrane density will influence the concentration-occupancy relationship. The consumption of the ligand would cause a shift to the right in the concentration response curve, with the size of the shift further increased by the reduction in the receptor concentration.
The time course of ligand binding to the receptor can also be modelled (Eq 1.5). For a soluble ligand, the ODE model shows that the time to equilibrium decreases with increasing ligand concentration, but is independent of the receptor level (S1A and S1B Fig). The stochastic nature of the interaction can be modelled using Gillespie's algorithm [19]. The degree of occupation varies with the association and dissociation constants and receptor number, with a greater variation in occupancy with fewer receptors expected due to the increase in stochasticity (S1C Fig). In summary, for a soluble monovalent ligand, the ODE modelling shows the relationship between ligand concentration, receptor occupancy and time, and how these vary in proportion to the association and dissociation constants and receptor number. The equilibrium receptor occupancy curve is approximately right-shifted when moving to the ligand-depletion model for a membrane-bound ligand.

ODE modelling of the interaction of a divalent ligand and a monomeric receptor
We have used ODEs to model the interaction of a divalent ligand with a monomeric receptor. In the equations below the divalent ligand is depicted as L 2 to reflect two epitopes and the receptor as R or RR to reflect a monomeric or dimeric receptor, respectively. There are two species of ligand-receptor combinations, L 2 R and L 2 RR. The introduction of a second epitope generates two equilibria: ð2:1Þ with equilibrium dissociation constants, K D1 and K D2 , respectively: https://doi.org/10.1371/journal.pcbi.1010708.g002 The two reactions (Eqs 2.1 and 2.2) have a shared intermediate species L 2 R acting as a product for Eq 2.1 and as a substrate for Eq 2.2. The relationship between the K D values of coupled reactions 1 and 2 is multiplicative [20] and introduces a new constant K D3 which is the product of K D1 and K D2 : ð2:5Þ The units of K D3 are M 2 , and so this constant is not equivalent to K D1 and K D2 which have units of M.
The binding of the first epitope in a soluble divalent receptor will influence the binding of the second epitope as a result of avidity, while dissociation from either receptor would return the dimeric receptor, L 2 RR, to L 2 R. A change in conformation of the divalent ligand could also influence the association and dissociation rate constants necessitating the introduction of an additional term, denoted by α, which would change the association and dissociation rate constants for the interconversion of L 2 R and L 2 RR: The ratio of α + and α − can be defined as a single variable, α, a cooperativity coefficient with a value of >1, equal to 1 and <1 denoting positive, neutral and negative cooperativity, respectively [21]: The value of α consequently influences the K D2 value, which is given by: The concentration-response relationship for the formation of L 2 R and L 2 RR with increasing concentrations of L 2 and its dependency on the association rate constant k 2 , while keeping the other rate constants the same, is shown in Fig 2B. For the purpose of modelling, the values of the dissociation rate constants k-1 and k-2 are the same as k 1 . The curves are generated by running the ODE solvers in Matlab to find the long-time solution. The curves can also be generated using an ODE model by solving a quadratic equation [22].
The amount of L 2 RR shows a bell-shaped relationship to the ligand concentration with the peak of L 2 RR occurring when the concentration of ligand is equal to K D1 and then declining as the ligand concentration increases, while the level of L 2 R continues to increase (Fig 2B). The magnitude of the peak in L 2 RR also increases with increasing values of k 2 . This is to be expected as this favours the formation of L 2 RR from L 2 R and reflects the influence of avidity. The modelling also shows that the ratio of L 2 R to L 2 RR is influenced by the receptor concentration, with a greater proportion of L 2 RR with increasing concentrations of the receptor. In other words, when the receptor concentration is low the probability of a divalent ligand finding a second receptor to form L 2 RR is reduced. The time to equilibrium also decreases as the concentration of ligand and receptors increases.
In summary, the introduction of a second epitope in the ligand leads to a further state, L 2 RR, and two new rate constants, k 2 and k -2 . This leads to a bell-shaped response curve for L 2 RR, which peaks when the concentration of ligand is equal to K D1 and with the size of the peak increasing in proportion to the forward rate constant, k 2 , while keeping all other rate constants the same. The amount of L 2 RR for a given ligand concentration also increases as the receptor number increases.

ODE modelling of the interaction of a tetravalent ligand and a monomeric receptor
Many of the ligands of platelet glycoprotein receptors have more than two epitopes such as extracellular matrix proteins, snake venom proteins, multisulfated sugars and diesel exhaust particles [23]. We have used a tetravalent ligand (depicted as L 4 ) as representative of a multivalent ligand and modelled this using ODEs. The interaction of the tetravalent ligand with a monomeric receptor will generate four ligand-receptor species: L 4 R, L 4 RR, L 4 RRR and L 4 RRRR, and therefore four equilibria defined by the order of attachment from unsaturated ligand to fully saturated ligand, with corresponding K D values for each reaction (Eqs 3.1-3.5 in S1 Text). Fig 2C shows the effect of the ligand and receptor concentration on the occupancy curve at equilibrium with neutral cooperativity (α = 1). The relationships between the concentration of the ligand and receptor occupancy, other than for occupancy of a single receptor, L 4 R, are bell-shaped with L 4 R predominating at higher ligand concentrations ( Fig 2C). In contrast, with increasing receptor concentrations L 4 RRRR will predominate for a fixed concentration of L 4 ( Fig 2C). Thus, as with a divalent ligand, the concentration of the ligand and the concentration of receptors are important in determining the stoichiometry of ligand-receptor binding as ligand valency increases.

Experimental validation of ODE models
Tetravalent but not divalent nanobodies activate GPVI. The ODE models define the relationship between receptor occupancy and ligand valency. To test the models, functional studies were undertaken using novel nanobody-based monovalent, divalent and tetravalent ligands to the platelet collagen receptor GPVI based on the assumption that a higher ligand valency would increase affinity (as a consequence of avidity) and signal strength. The divalent and tetravalent ligands were generated by crosslinking of a high affinity nanobody, Nb2, that binds to the first immunoglobulin (Ig) domain in GPVI [24]. A flexible 15 amino acid linker (GGGGS) 3 was used to generate a divalent ligand, Nb2-2, and a tetravalent ligand, Nb2-4, as shown in S2A Fig Neither Nb2 nor Nb2-2 (150 nM) induce platelet aggregation when incubated for up to 60 min ( Fig 3A[i]). On the other hand, Nb2 and Nb2-2 block aggregation to collagen (3 μg/mL) with IC50 values of 30.7±1.9 nM and of 4.7±0.5 nM, respectively ( Fig 3A[ii-iv]). Nb2-2 also exhibited a 6-fold greater affinity (K D = 0.100 + 0.003 nM) than Nb2 (K D = 0.58 + 0.06 nM) when measured by surface plasmon resonance (SPR; S3 Fig and [24]). The greater potency and affinity of Nb2-2 relative to Nb-2 shows that both epitopes can bind to GPVI on platelets and on immobilised protein, respectively, with the increase in affinity being due to avidity.
Nb2-4 stimulated rapid aggregation of platelets with a maximal response at a low nanomolar concentration suggesting that it is able to crosslink at least 3 and possibly 4 GPVI receptors (Figs 3B[i-ii] and S2C [i,ii]). The increase in crosslinking accounts for its ability to induce sufficient clustering of GPVI to drive platelet aggregation. As expected, aggregation to Nb2-4 was blocked by monovalent Nb2 and by inhibitors of Src and Syk tyrosine kinases (Fig 3[iii-iv]) confirming that it was mediated through clustering of GPVI. The highest concentration of Nb2-4 tested (100 nM) induced full platelet aggregation ( Fig 3B[ii]) with no evidence of a bellshaped concentration response relationship as seen in the ODE model. The affinity of Nb2-4 for GPVI measured by SPR was similar to that of Nb2-2 (S3 Fig). These results demonstrate that a tetravalent but not a divalent ligand is able to activate platelets through GPVI and that a divalent ligand has a greater potency than a monovalent ligand for inhibition of platelet activation by collagen.

Agent-based modelling
Agent-based modelling to monitor ligand binding and cluster formation. We have developed an agent-based model in NetLogo [25] to complement the modelling of ligand- For the agent-based model, the following assumptions apply: (i) the crosslinker is only able to bind to a phosphorylated receptor (ii) the mobility of monomeric and dimeric receptors is the same (iii) all reactions are reversible (iv) the speed of a cluster (> 3 receptors) is inversely proportional to its size (v) the binding of the cytosolic crosslinker requires receptor phosphorylation and protects against dephosphorylation. Receptor clustering can be achieved by a combination of (i) receptor-receptor homo-dimerisation (ii) binding of a cytosolic crosslinking protein (iii) a divalent or multivalent ligand. Together, these can lead to formation of receptor trimers, tetramers and higher order oligomers. The assumption that the binding of the cytosolic crosslinker requires receptor phosphorylation and protects against dephosphorylation reflects the situation in a cell, as shown using phosphorylated peptides [26]. The reduction in diffusion of large clusters was introduced as they diffuse at a slower rate in the membrane [27], but when tested in a limited number of scenarios, removal of this parameter did not have a significant effect on protein phosphorylation or receptor clustering. The default values for on-rates and off-rates are shown in ticks (arbitrary unit of time), with one tick corresponding to the 'go' function within each run. These represent the probability of a reaction to occur out of 100 runs. These are arbitrary values selected for illustrative purposes. They have been chosen to give a low basal level of dimerisation and phosphorylation to correspond to a resting platelet. Two sets of values have been chosen to reflect a moderate and high affinity crosslinker in regulating cluster formation; the latter takes into account that the affinity of a 'moderate' affinity crosslinker can be increased by its co-localisation in the membrane as is the case for the majority of SH2-domain containing proteins. A 'low' value was selected for the affinity of the ligand to reflect the affinity of individual epitopes in endogenous ligands for platelet glycoprotein receptors, noting that this can be increased by avidity. The selection of a relatively low value enables illustration of how variables can combine to drive clustering. The rate of movement of receptors is constant at 1 patch per tick patch is an arbitrary unit of area (there are 441 patches per box). The default values give rise to to a basal value of 10% receptor dimerisation and 10% receptor phosphorylation in the absence of crosslinker. The crosslinker can only bind to a phosphorylated receptor: the probability of each epitope of the crosslinker to bind to a phosphorylated receptor is shown as an average frequency per 100 ticks. Assuming that both epitopes are occupied, the 'moderate' affinity crosslinker will remain bound to two receptors for 16 out of 100 ticks on average. The 'high' affinity crosslinker will remain bound to the receptor for 81 out of 100 ticks, on average. The probability of each epitope of the divalent ligand to attach to a receptor upon collision was set to 20 out of 100 ticks and to dissociate to 20 out of 100 ticks. Assuming that both epitopes of the divalent ligand are bound at the same time, a fully occupied divalent ligand will remain bound to both receptors for an average 4 out of 100 ticks. https://doi.org/10.1371/journal.pcbi.1010708.g004 We have systematically studied each of these parameters and show their effect on receptor dimerisation, receptor clustering and receptor phosphorylation.
Agent-based modelling of receptor dimerisation in the absence of a ligand or crosslinker. We first used ABM to model dimerisation of a receptor that undergoes reversible dimerisation in the absence of ligand or a crosslinker (Eq 4.1): ð4:1Þ The model shows that the degree of dimerisation is influenced by the receptor number and the rates of association and dissociation. As the receptor number increases, the number of collisions between receptors increases leading to an increase in the proportion of dimers increases (Fig 5A[i-iii] and S1 Video). Increasing the rate of association or decreasing the rate of dissociation also increases the number of receptor dimers (Fig 5B and 5C Agent-based modelling of receptor dimerisation in the presence of a divalent crosslinker. Introducing a divalent crosslinker that undergoes reversible binding to a phosphorylated receptor provides an alternative mechanism of clustering. In the agent-based model, the binding of the crosslinker prevents receptor dephosphorylation leading to a net increase in phosphorylation. This mimics the effect of binding a SH2 domain to a phosphorylated tyrosine in a cell. The effect of the crosslinker was modelled with on-and off-rates of receptor phosphorylation set to achieve 10% basal phosphorylation at the initiation of the run. This level of phosphorylation was selected to reflect the low level of receptor phosphorylation in a resting cell. The model was run in the presence of a 'moderate' or a 'high' affinity crosslinker as defined using the parameters in whereas the presence of a high affinity crosslinker led to an increase in phosphorylation to over 50% and a doubling in the number of receptor dimers. There was no formation of higher order clusters, as expected, as dimerisation of the receptor monomers can only occur through binding of the crosslinker (Fig 6A[i-iii] and S2 Video).
The model was then run for a receptor that is able to undergo reversible dimerisation. The basal value of dimerisation was set at 10% in the absence of crosslinker (Fig 4). The introduction of a moderate affinity crosslinker resulted in a small increase in the number of receptor dimers and doubling of receptor phosphorylation at the highest crosslinker concentration, but there was negligible formation of higher order clusters (S4B[i-iii] Fig and S3 Video). The presence of a high affinity crosslinker led to an increase in the number of receptor dimers and the degree of receptor phosphorylation (both values reaching over 50%) and a small increase in the number of receptor trimers and tetramers (Fig 6B[i-iii]). The increase in the number of dimers, trimers and tetramers reflects the crosslinking of receptor monomers and dimers to themselves and to each other.
The model was extended to investigate the effect of increasing the basal level of receptor phosphorylation and receptor dimerisation. For these studies, the number of crosslinkers was set to the default parameter of 100 per box. The 1: 1 ratio of receptors to crosslinker is similar to the ratio of CLEC-2 and Syk in human platelets measured by quantitative proteomics (S1 Table). For comparison, S1 Table also shows the levels of these proteins in mouse platelets, as well as the levels of FcγRIIA, GPVI and PEAR1 in human and mouse platelets.
Increasing the rate of receptor phosphorylation had a negligible effect on the number of dimers and trimers in the presence of a moderate affinity crosslinker (S4C[i-iii] Fig). In contrast, in the presence of a high affinity crosslinker over 50% of receptors were present as dimer, trimers, tetramers and higher order oligomers (Fig 6C[i- These results show that receptor dimerisation, phosphorylation and crosslinker affinity regulate the degree of clustering, and that in combination they can lead to significant levels of clustering.
Agent-based modelling of receptor dimerisation in the presence of a divalent ligand. The readouts of the agent-based models are receptor clustering and phosphorylation, neither of which are altered in the presence of a monovalent ligand. For this reason, the models would not change in the presence of a monovalent ligand.
We have therefore used agent-based model to investigate the effect of a divalent ligand mimicking for example an antibody or dimerised nanobody. As this is similar to modelling the effect of a divalent crosslinker with a high basal level of receptor phosphorylation, as shown in Fig 6D, the modelling has only been carried out in combination with a divalent crosslinker. Low basal values of receptor dimerisation and phosphorylation were used, along with a moderate affinity crosslinker, in anticipation that the combination of these parameters may lead to synergy.
We first modelled the effect of a divalent ligand with a monomeric receptor that is unable to dimerise in the presence of a moderate affinity crosslinker and a basal level of receptor phosphorylation of 10%. This led to an approximate doubling in receptor phosphorylation and receptor dimerisation and formation of a small number of trimers and tetramers at the highest number of ligands (Fig 7A[i-iii]). In the presence of a low basal level of receptor dimerisation and absence of crosslinker, the divalent ligand led to a similar increase in receptor phosphorylation along with the formation of receptor trimers, tetramers and higher order oligomers even in the presence of a low ligand number (Fig 7B[i-iii]). Strikingly, the further addition of the moderate affinity crosslinker led to the formation of large receptor clusters which in some runs resulted in the formation of a single cluster (Fig 7C[i-iii] and 7[iv] and S6 Video).
This analysis shows that a divalent ligand can induce formation of a large receptor clusters in the presence of a low basal of receptor dimerisation and moderate affinity crosslinker. The relatively low nature of these values would help to minimise constitutive signalling in a cell whilst enabling rapid amplification upon ligand engagement.

Experimental validation of the agent-based models
The effect of a Syk inhibitor on clustering of CLEC-2 and GPVI. The agent-based model in Fig 7C shows how the combination of a divalent ligand, receptor dimerisation, and crosslinker can drive receptor clustering. This model is applicable to CLEC-2 and PEAR1 whose phosphorylated tails are crosslinked by the tandem SH2 domains in Syk and PI 3-kinase, respectively. This is in contrast to GPVI where the tandem SH2 domains in Syk are believed to bind to two phosphorylated tyrosines in the FcR γ-chain ITAM.
We used fluorescence correlation spectroscopy (FCS) to investigate a role for Syk in the clustering of CLEC-2 and GPVI. FCS measures the fluorescence fluctuations brought about by the diffusion of fluorescently-tagged proteins through a small volume (~0.2 fl) in the membrane, generating information on the rate of diffusion and molecular brightness. The change in diffusion is relatively insensitive to cluster formation as an 8-fold change in mass results in a 1.6-fold change in diffusion coefficient [27]. In contrast, when used with photon counting histogram (PCH) analysis, FCS can be used to differentiate monomers, dimers and higher order oligomers [28].
To investigate the role of Syk in the clustering of CLEC-2, we expressed a N-terminal eGFP-tagged CLEC-2 in HEK293T cells and stimulated the cells with the snake venom toxin rhodocytin and the monoclonal antibody (mAb) to CLEC-2, AYP1. FCS measurements were carried out in the absence and presence of the Syk inhibitor, PRT-060318, which blocks phosphorylation of the hemITAM in CLEC-2 [29]. The membrane localisation of the eGFPtagged CLEC-2 was verified by confocal microscopy (Fig 8B[i]).
Rhodocytin and mAb AYP1 induce clustering of CLEC-2 on the surface, as shown by confocal microscopy, which are blocked in the presence of PRT-060318 (Fig 8B[i]). In line with this, both ligands increase the molecular brightness of CLEC-2 as measured by FCS and this is also blocked in the presence of PRT-060318 (Fig 8B[ii-iii]) confirming a critical role of Syk in clustering of CLEC-2.
In contrast, PRT-060318 had no effect on clustering of C-terminal eGFP-tagged GPVI measured by confocal microscopy and FCS upon activation by a collagen-related peptide (CRP) or by the divalent and tetravalent nanobodies, Nb2-2 and Nb2-4 (S5A and S5B[i,ii] Fig). This is consistent with the tandem SH2 domains in Syk binding to the phosphorylated ITAM and not crosslinking adjacent receptors.
Synergistic action of a blocking F(ab) 2 and Syk inhibitor in blocking platelet activation. The demonstration that clustering is regulated by binding of a divalent ligand and by a moderate affinity crosslinker suggests that blocking both pathways may have a synergistic effect. To investigate this, we tested the combination of a CLEC-2-blocking antibody fragment, AYP1 F(ab) 2 (2.2 nM), and the Syk inhibitor, PRT-060318 (2.5 nM), at concentrations that have a minor/negligible effect on the response to rhodocytin. In combination, they caused a marked inhibition of the response to rhodocytin (Fig 9A[i,ii]).
To investigate if the inhibition could also reflect loss of signalling by Syk, we performed a similar study with the GPVI blocking nanobody, Nb2-2, and PRT-060318 (2.5 nM) against CRP (1 μg/ml). Nb2-2 has a similar subnanomolar affinity for GPVI to that of AYP1 (Fab 2 ) for CLEC-2, and was used at the same concentration (2.2 nM) except in one donor whose platelets exhibited a much greater sensitivity to inhibition. The combination of the two inhibitors caused a marked inhibition of response to CRP (1 μg/ml), although in some experiments this resulted in a delay rather than abolition of aggregation, with a maximal response observed at later times (Fig 9B[i,ii]). This was not seen in the studies with rhodocytin.

Discussion
The binding of multivalent ligands to platelet glycoprotein receptors such as CLEC-2, FcγRIIA, GPVI and PEAR1 leads to dimerisation and higher order clustering which drives formation of a critical density of tyrosine-based signalling motifs in the inner leaflet of the membrane for phosphorylation by active, membrane-associated kinases. This results in the binding of SH2 domains thereby overcoming the effect of active membrane phosphatases. In the present study, we have used ODE and ABM to model the interaction of these receptors with their ligands on the assumption that the interaction is driven by the law of mass action and that additional clustering can be achieved through receptor homo-dimerisation and a tandem SH2 domain protein.
In the ODE models, we have assumed that each epitope in a ligand binds to one receptor, that the concentration of the ligand is in vast excess of the receptor (unless stated), and that receptors are monomers. The ODE models show a bell-shaped concentration-response relationship for binding of divalent and tetravalent ligands (the latter representing multivalent ligands). In the ODE models, the divalent ligands binds to mixture of monomers dimers, and the tetravalent ligands to mixture of monomers, dimers, trimers and tetramers. Because of the effect of avidity, divalent and tetravalent ligands have a greater affinity for the receptor, as illustrated by SPR measurements for the divalent and tetravalent nanobody ligands to GPVI relative to the monovalent backbone. Further, confocal microscopy and FCS studies show that the two multivalent ligands but not the monovalent backbone cluster GPVI in a transfected cell line, thereby showing that they can crosslink GPVI on a cell surface. Strikingly, however, Nb2 and Nb2-2 are antagonists of collagen-induced activation of GPVI in human platelets, whereas tetravalent Nb2-4 induces powerful platelet aggregation. This suggests that four and possibly as low as three epitopes in a ligand are required for sufficient clustering of GPVI to drive platelet activation.
The concentration-response curve to Nb2-4 did not exhibit a bell-shaped relationship at higher concentrations in contrast to the ODE models. This is likely to be due to a combinaton of factors including (i) the stochastic nature of the ligand-receptor interaction, (ii) full receptor occupancy is not required for maximal aggregation by GPVI as shown by stimulation of platelets by CRP [30,31] and (iii) GPVI is expressed as a mixture of monomers and dimers as a consequence of homo-dimerisation and co-location on the surface [32]. The presence of dimers will lead to larger receptor complexes through crosslinking of monomers-to-dimers and cells (field of view = 52 x 52 μm) (scale bar = 5 μm). (Bii) Box plots showing the molecular brightness (cpm s -1 ) of CLEC-2-eGFP, centre lines represent median, box limits indicate the 25 th and 75 th percentiles, and whiskers extend to minimum and maximum points. Concentration of rhodocytin (30 nM), mAb AYP1 (6.6 nM), PRT (10 μM) were used. Significance was measured with Kruskal-Wallis with Dunn's post-hoc where P � 0.05. In (ii) # = significance compared to CLEC-2 alone. FCS measurements were taken in 35-48 cells; n = 3-5.
https://doi.org/10.1371/journal.pcbi.1010708.g008 dimers-to-dimers. While the presence of dimers can be modelled using ODE, this rapidly becomes complex and this does not cover co-association by chance. On the other hand, spatial consideration can be readily modelled by ABM, as well as other factors including receptor dimerisation and a cytosolic crosslinker (mimicking a tandem SH2 domain protein).
The agent-based model developed in this study shows how the presence of a divalent ligand, receptor homo-dimerisation and cytosolic crosslinker (and receptor phosphorylation) can work in combination to generate large receptor clusters. The importance of the crosslinker in CLEC-2 signalling was verified experimentally using confocal microscopy and FCS using an inhibitor of Syk which blocks phosphorylation of the CLEC-2 hemITAM [29]. Syk has a two SH2 domains which crosslink two phosphorylated CLEC-2 receptors. The agent-based model also shows a correlation between receptor density and receptor dimerisation, which may explain the increase in NFAT signalling that is observed with increasing GPVI concentrations in a transfected DT40 B cell line in the absence of a ligand, as this is predicted to increase receptor dimerisation [33]. This indicates that dimerisation of GPVI is sufficient to drive signalling. In contrast, the divalent nanobody against GPVI, Nb2-2, was an antagonist in human platelets. Similarly, an antibody to FcγRIIA, mAb IV.3, is an antagonist in human platelets [34]. The difference between the transfected cell line model and platelets may be related to the the level of GPVI expression and/or the levels of active kinases and phosphatases in the membrane, as well as differences in the sensitivity of the two assays. In platelets, there is a high level of basal tyrsoine phosphatase activity, as shown by the rapid increase in tyrosine phosphorylation in the presence of the phosphatase inhibitor peroxovanadate [35], and this helps to minimise activation of ITAM and YxxM receptors in the absence of a strong activating ligand. The importance of receptor diffusion in supporting receptor dimerisation and higher order clustering accounts is likely to account for the delay in onset of platelet activation by multivalent ligands of CLEC-2 [36], GPVI [37] and FcγRIIA [38,39].
The ODE model was also applied to the interaction of a membrane or immobilised ligand. In both cases, the binding of the ligand to a receptor leads to a reduction in its concentration and thus a shift in the concentration-response curve to the right, although in reality this is more than countered by avidity. In the case of an immobilised multivalent ligand, clustering is achieved through the mobility of the receptor on the platelet surface as illustrated for the binding of platelet GPVI to collagen, as shown using single molecule super-resolution microscopy [40]. On the other hand, the clustering of CLEC-2 would require either a degree of prior clustering of podoplanin in the membrane or the formation of clusters as a result of binding and diffusion. The clustering of podoplanin could be achieved through binding of its cytosolic tail to the cytoskeletal via the ezrin, radixin and moesin (ERM) family of binding proteins [5] or in the case of a high level of expression by diffusion. The clustering and phosphorylation of CLEC-2 receptors to the podoplanin-rich areas would be reinforced by the binding of Syk to the CLEC-2 phosphorylated tail, thereby bringing in additional CLEC-2 receptors and further clustering of podoplanin. This has been demonstrated in lipid bilayers, with clustering of podoplanin by platelet CLEC-2 being dependent on Src and Syk tyrosine kinases [41]. The combined effect of avidity, diffusion and crosslinking by Syk counters the need for a high affinity interaction between podoplanin and CLEC-2 in mediating activation. In line with this, we have estimated the concentration of CLEC-2 in the membrane using to be below the reported affinity constant (K D ) of 24 μM for binding to podoplanin [42] (see S1 Text). The net effect of affinity and avidity necessitates that competitive inhibitors would require a high affinity and a slow rate of dissociation for effective and sustained antagonism.
There are several factors that limit clustering of receptors under basal conditions and therefore unwanted activation. For platelet glycoprotein receptors, this includes their predominantly monomeric expression as illustrated for GPVI [32] and a low constitutive level of receptor phosphorylation as a result of a high tyrosine phosphatase activity at the membrane [6,35,43]. Receptors may also be prevented from dimerising by low affinity interactions with other membrane proteins, notably tetraspanins which have little biological activity on their own but have been shown to associate with platelet surface receptors [44]. It follows that a change in the distribution or activity of proteins that associate with or are located close to GPVI, such as tetraspanins and tyrosine phosphatases, has the potential to drive activation of glycoprotein receptors in the absence of ligand engagement.
There are a number of limitations to be considered in the models that have been developed in this study. In the ODE models, the single transmembrane receptors are considered to be monomeric, but in a cell there will be a mixture of dimers and higher order oligomers as discussed above for GPVI. The ODE and agent-based models are based on the presence of a single receptor for each ligands, but several ligands including collagen and fibrin bind to more than one receptor on the platelet surface. The modelling of a ligand with two different receptors is beyond the scope of this study.
In summary, the models developed in this study highlight parameters that influence occupancy and clustering of weakly dimerising glycoprotein receptors, including ligand valency, receptor level and a cytosolic crosslinker. Predictions of the models have been tested through generation of novel divalent and tetravalent ligands, and combined use of inhibitors of ligand binding and Syk in functional studies. The results show that a tetravalent but not a divalent ligand induces powerful activation of human platelets demonstrating that 3 and possibly 4 ligands are required to induce activation. The absence of a bell-shaped dose response-relationship for the tetravalent ligand, in contrast to the ODE equlilbrium binding relationships, is like to reflect the existence of of pre-existing receptor dimers and other factors as discussed above. Fig 10 illustrates how receptor dimerisation and a cytosolic crosslinker combine to drive Fig 10. A model depicting how a combination of a weakly dimerising receptor and a crosslinking agent can drive cluster formation. The key is shown on the upper right: (A) illustrates crosslinking of two phosphorylated receptors in such a way that growth can only occur through dissociation or addition of a divalent ligands (B) ligand-induced dimerisation of two crosslinked dimers (C) cluster growth by a dimerising receptor and crosslinker. This can be further increased in the presence of a divalent ligand.
https://doi.org/10.1371/journal.pcbi.1010708.g010 clustering, and how this can be increased in the presence of a divalent ligand. The synergy between low concentrations of a receptor blocker and a kinase inhibitor provides a strategy to overcome the effect of avidity and to lower the off-target effects of inhibitors.
There have been relatively few studies modelling activation of the platelet glycoprotein receptors that are the focus of this study. Martyanov and colleagues [45] have developed a stochastic multicompartmental computational model of CLEC-2 signalling which identifies clustering and signalosome formation as critical rate-limiting steps in platelet activation. The loss of downstream signalling accounts for the synergy between threshold concentrations of a GPVI-blocking nanobody and Syk inhibitor, even though Syk does not play a role in clustering of GPVI. The ability of Syk inhibition to block of CLEC-2 accounts provides an additional basis for the synergy and may explain why the inhibition does not recover at later times. Griffie et al [46] have modelled the clustering of receptors using ABM using a pre-set rule of "desire for clustering". The application of this rule in the model leads to discrete clusters of receptors of uniform size (assuming a starting random distribution of the receptor). This study together with the models developed in thie manuscript illustrate how ABM can lead to a new understanding of the mechanisms of receptor clustering based on simple rules which can be tested in wet-lab experiments.

Ethics statement
Healthy human volunteers for blood donation gave written consent under the Programme of Work ERN_11-0175 "The regulation of activation of platelets" ERN_11-0175P which was approved by the Science, Technology, Engineering and Mathematics Ethical Review Committee of the University of Birmingham.

Modelling methods
Deterministic and stochastic modelling using ordinary differential equations (ODEs). Systems of equations based on the law of mass action were developed to model the binding of monovalent, divalent and tetravalent soluble ligands, and monovalent membrane-bounds ligands to monomeric receptors. The numerical solutions for deterministic simulations were calculated using MATLAB's in-built ODE solver ode15s, which was chosen for its advantage of being able to compute solutions for stiff differential equations (version 2019a). To compute stochastic simulations, the approach implemented was based on Gillespie's algorithm [19]. Random number generators were used to introduce a degree of randomness for the time when a reaction would take place while considering the probability of each reaction occurring based on the rate constants, and the probability of the next reaction to occur was weighted by that of the previous reaction in addition to its respective rate constant. The outcome of x1000 iterations was combined and plotted against time.
The scripts for developing ODE models can be found at: https://github.com/zeemaqsood/ TAPAS_ESR9_Modelling_Project

Agent-based modelling
ABM for monovalent, divalent and tetravalent ligands was performed using NetLogo [25] and used to study the interaction of such ligands with a monomeric receptor. The ABM models comprise of agents (which are termed 'turtles' in NetLogo) of multiple types, or breeds, and their interactions within the modelling environment based on defined rules. There are three main breeds: ligands, receptors and cross-linking proteins in the cytosol. Each breed possesses a unique set of attributes that can be changed and are responsible for conversion from one state to another as explained in the S1 Text and S6 Fig. Each run of 3,000 ticks was repeated 30 times, the results are shown as mean + s.d. The scripts for ABM model used can be found at: https://github.com/zeemaqsood/TAPAS_ESR9_Modelling_Project

Experimental methods
Materials. Rhodocytin was purified from the venom of Calloselasma rhodostoma as described [47]. Atto-488 dye was obtained from Merck Life Science UK Limited (Dorset, UK). The Syk inhibitor PRT-060318 was purchased from Caltag Medsystems (Buckingham, UK). CRP was purchased from Cambcol Laboratories (Cambridge, UK). Other chemicals/reagents were obtained from Merck Life Science UK Limited or Polysciences (Pennsylvania, USA).
Antibodies and nanobodies. Nanobodies were raised against GPVI extracellular region through VIB Nanobody core (VIB Nanobody Service Facility, Brussels, Belgium). The most potent of these, Nb2, was expressed and purified as previously described [24]. Divalent (Nb2-2) and tetravalent (Nb2-4) versions of Nb2 were cloned by VIB Nanobody Core to contain a flexible (GGGGS) 3 linker [48,49] between each individual Nb2 domain, see S2A Fig for domain architecture. Expression and purification of Nb2-2 and Nb2-4 was as described for monovalent Nb2: E.coli WK6 cells were used to express the Nb construct which contained an N-terminal PelB leader sequence to secrete protein into the periplasmic space. The constructs also contained a C-terminal His-6 tag which was used for purification by nickel affinity chromatography. The mouse mAb AYP1 was raised against the C-type lectin-like domain of human CLEC-2 and Fab and F(ab) 2 fragments generated as described [50].
Platelet preparation and light transmission aggregometry. Blood was obtained from healthy volunteers who had denied taking antiplatelet medication in the previous 10 days using 4% sodium citrate.
Human CLEC-2-eGFP and GPVI-eGFP construct generation. cDNA for human CLEC1b (the gene encoding CLEC-2) with a linker coding for DRNLPPLAPL was synthesised (GeneArt, Invitrogen). The N-terminally eGFP tagged CLEC-2 construct was made by subcloning Clec1b cDNA in the pEGFP(A206K)-C1 expression vector for expression of CLEC-2-eGFP where the N-terminal eGFP has a A206K mutation to prevent dimerisation of the fluorescent protein. Generation of the GPVI-eGFP construct has been described previously [32].
Fluorescence correlation spectroscopy (FCS) and analysis. HEK293T cells were seeded in phenol red-free DMEM at a density of 3x10 4 cells/25mm coverslip (Marienfeld, high precision, thickness No. 1.5H (0.170 mm ± 0.005 mm). The next day, cells were transiently transfected using polyethylenimine (PEI Max MW 40,000, Polysciences) (PEI:DNA ratio = 3:1, 3 μl:1 μg) where 100 ng CLEC-2-eGFP or 100 ng GPVI-eGFP + 100 ng FcRγ-chain DNA was used to achieve optimal receptor density. Cells were left growing a further 24 h at 37˚C/5% CO 2 . FCS measurements were made using a Zeiss LSM-880 confocal microscope equipped with gallium arsenide phosphide photon detectors (GaAsP) (Carl Zeiss, Jena, Germany) and prior to each experiment the microscope was aligned and calibrated using Atto-488 dye as described previously [32]. On the day of the experiment, the FCS confocal volume was determined by measurement of the axial and lateral radii. For FCS experiments with ligands, rhodocytin (30 nM), AYP1 mAb (6.6 nM), CRP (10 μg/ml), Nb2-2 (10 nM) and Nb2-4 (10 nM) were made in phenol red-free DMEM and added to the cells and imaged immediately. For PRT samples, PRT-060318 (10 μM) was made in phenol red-free DMEM and added to cells for 45 min prior to ligand addition and imaging. DMSO only (0.05 and 0.1%) samples were imaged as a control. FCS data were analysed by photon counting histogram (PCH) analysis to determine molecular brightness using Zen 2012 (black edition) software (Carl Zeiss, Jena, Germany) as described previously [32].
Surface plasmon resonance. Surface plasmon resonance experiments were performed using a Biacore T200 instrument (Cytiva). GPVI (extracellular region) was immobilised directly on the CM5 chip using amine-coupling to the carboxylmethylated dextran-coated surface. Reference surfaces were blocked using 1M ethanolamine pH 8. All sensograms shown are double reference subtracted and at least two replicates were injected per cycle as well as experimental replicates of n = 3. Experiments were performed at 25˚C with a flow rate of 30 μl/min in HBS-EP running buffer (0.01 M HEPES pH 7.4, 0.15 M NaCl, 3 mM EDTA, 0.005% v/v surfactant P20). Multi-cycle kinetic assays were used with at least 5 concentration points between 0.1x and 10x the K D . Each concentration of analyte (GPVI nanobody divalent/tetravalent) was run as follows, 120 sec injection, 900 sec dissociation, 30 sec regeneration with 10 mM glycine pH 1.5 followed by a 300 sec stabilisation period. Kinetic analysis was performed using the Biacore T200 Evaluation software using a global fitting to a 1:1 binding model.

Statistical analysis
Results are shown as mean ± s.d. unless otherwise stated and the number of independent experiments is described in Figure legends. Dose response curves were fitted using doseresponse fit in Origin (version 8.6, OriginLab, Northampton, MA). For FCS, data sets were first tested for normality using the Shapiro-Wilks test. FCS data were analysed by Kruskal-Wallis with Dunn's post-hoc test using PRISM v9.2.0 (GraphPad, San Diego, CA). For the effect of inhibitors on aggregation, all treatment groups were first compared with each other using one-way ANOVA analysis, followed by multiple comparison tests. Significance was set at P � 0.05. There are six reactions occurring on three species: a divalent ligand (L 2 ), receptor (R) and a crosslinker (S). A ligand-receptor complex formation reaction is coloured red, receptor phosphorylation (either basal or post ligand attachment to a receptor) is coloured blue, receptor-crosslinker complex formation is coloured green (either to an unbound or partially bound crosslinker), and receptor-receptor dimerisation reaction is coloured yellow. (TIF) S1 Table. Protein expression levels in platelets. The expression levels of monomeric receptors and Syk. The levels are taken from quantitative proteomic studies in human [51] and mouse platelets [52]. (TIF) S2 Table. Turtle breeds and their attributes. All turtles are differentiated into three main breeds namely ligands, receptors and cytosolic cross-linkers. Each breed possesses a unique set of attributes that can be changed. The attributes are responsible for conversion of the state of each agent. The main attributes for each breed are listed in the Table below