Ligand Binding Dynamics for Pre-dimerised G Protein-Coupled Receptor Homodimers: Linear Models and Analytical Solutions

Evidence suggests that many G protein-coupled receptors (GPCRs) are bound together forming dimers. The implications of dimerisation for cellular signalling outcomes, and ultimately drug discovery and therapeutics, remain unclear. Consideration of ligand binding and signalling via receptor dimers is therefore required as an addition to classical receptor theory, which is largely built on assumptions of monomeric receptors. A key factor in developing theoretical models of dimer signalling is cooperativity across the dimer, whereby binding of a ligand to one protomer affects the binding of a ligand to the other protomer. Here, we present and analyse linear models for one-ligand and two-ligand binding dynamics at homodimerised receptors, as an essential building block in the development of dimerised receptor theory. For systems at equilibrium, we compute analytical solutions for total bound labelled ligand and derive conditions on the cooperativity factors under which multiphasic log dose–response curves are expected. This could help explain data extracted from pharmacological experiments that do not fit to the standard Hill curves that are often used in this type of analysis. For the time-dependent problems, we also obtain analytical solutions. For the single-ligand case, the construction of the analytical solution is straightforward; it is bi-exponential in time, sharing a similar structure to the well-known monomeric competition dynamics of Motulsky–Mahan. We suggest that this model is therefore practically usable by the pharmacologist towards developing insights into the potential dynamics and consequences of dimerised receptors.


Introduction
Mathematical models have long played an important role in analytical pharmacology, which has its roots in receptor theory (Kenakin 2009;Kenakin and Williams 2014;Kenakin and Christopoulos 2011;Kenakin 2004;Colquhoun 2006). Much of this theory rests upon assumptions of equilibrium and ligands binding to monomeric receptors. It is now widely accepted that consideration of binding and signalling dynamics is an important factor in the drug discovery process (Schuetz et al. 2017). Furthermore, there is widespread acknowledgement that many receptors may exist (constitutively) as dimers or higher oligomers (Milligan 2004(Milligan , 2006(Milligan , 2013Milligan and Smith 2007). Therefore, a theoretical foundation for the study of signalling via dimerised receptors is required in order to classify, quantify and simulate ligand-receptor interactions and their signalling outcomes, for both equilibrium and dynamic conditions. G protein-coupled receptors (GPCRs) are membrane receptors of particular therapeutic importance, since they represent targets for up to approximately 50% of current drugs ). These receptors were long thought to exist solely as monomers, and theoretical models for GPCR signalling have largely been built around this assumption, in line with classical receptor theory for both equilibrium and dynamic models (Woodroffe et al. 2009(Woodroffe et al. , 2010Bridge 2009;Bridge et al. 2010). However, it has been demonstrated and is now widely recognised that GPCRs may exist and function as dimers and higher-order oligomers (Smith and Milligan 2010;Milligan 2004Milligan , 2006Bai 2004). The extent to which a ligand bound to one site of an oligomeric receptor complex affects ligand binding characteristics to other sites is termed cooperativity, and the effect of cross-dimer (or cross-oligomer) cooperativity on binding and signalling dynamic responses is yet to be fully elucidated theoretically.
In pursuit of potentially developing drugs which exploit ligand binding and activation crosstalk at dimeric receptors, many questions remain. For example, what are the potential advantages (to a cell or tissue) to express dimeric rather than monomeric receptors? "Expanded pharmacology" possibilities (in terms of binding and activation) may increase the diversity of potential signalling responses (via G protein coupling and signalling) towards "fine-tuning" cell responses (Smith and Milligan 2010). Identifying homodimers and heterodimers and understanding their role in signalling is a key step towards exploiting oligomerisation in drug discovery (Smith and Milligan 2010;Milligan 2004). The prevalence of GPCR dimers is a subject of debate, but their capacity for functional outcomes beyond those offered by monomers gives clear potential for the development of novel therapeutics (Milligan 2013). Indeed, with respect to classical receptor theory, dimerised receptors may represent "novel receptors" , i.e. receptors whose binding and signalling response profiles are more complex and offer more targetable therapeutic potential than their constituent protomers alone. Much needs to be established in terms of differential pharmacology and functional effects before dimers become tractable drug targets (Milligan 2009); development of the theory of dimerised GPCRs is therefore needed (Bai 2004), and we begin here by presenting mathematical models of ligand binding dynamics.
Mathematical modelling of binding to dimerised GPCRs has largely focused on equilibrium binding models Chidiac et al. 1997;Durroux 2005;Franco et al. 2006Franco et al. , 2005Franco et al. , 2007Ferré et al. 2014), but the dynamics of ligand binding at such dimers have not been widely modelled. Typically, the equilibrium models yield analytical solutions for total bound labelled ligand, derived algebraically using mass action and receptor conservation considerations. Corresponding log dose-response (logDR) curves for bound ligand may be multiphasic, exhibiting multiple inflections (Ferré et al. 2014;Casadó et al. 2009;Rovira et al. 2009;Durroux 2005;Chidiac et al. 1997). For models of dimerised receptor binding, this departure from monophasic logDR curves (typical of monomeric receptor binding) may be quantified by a dimer cooperativity index (Giraldo 2008;Casadó et al. 2007;Franco et al. 2008), which relates to the apparent binding cooperativity in the Hill function sense. Within a Hill function analysis, however, there is the possibility that information is missed due to the inability of Hill coefficients to distinguish interaction mechanisms (Prinz 2010). These works appear to be the state of the art in practical GPCR models, while more mathematically abstract approaches are taken elsewhere: an algebra of dimerisation is presented in Woolf and Linderman (2004), and generalised multi-site binding models are analysed in Juška (2008).
Dynamic models of binding and signalling for dimerised receptors are less common than equilibrium models. Spatial models of the dimerisation process and subsequent signalling are developed in Mayawala et al. (2006), while dynamics of receptor and transducer protein dimerisation are studied using an ordinary differential equation (ODE) model in Vera et al. (2008), wherein it is suggested that homodimerisation may serve to regulate signalling over multiple timescales. Ligand-induced dimerisation of VEGF receptors is modelled using ODEs in Mac Gabhann and Popel (2007). In order to lay the foundations for dynamic modelling of signalling via dimerised GPCRs, we refer to a more recent study of GPCR ligand binding dynamics (May et al. 2011). Therein, linear ODE models of single-ligand and competition dynamics at pre-dimerised homodimeric receptors are presented, with a brief model analysis and numerical data fitting, without presenting analytical solutions to the ODE systems. In the current work, we describe a simplified formulation of the May et al. models and derive their analytical solutions. For the single-ligand model, the solution structure (bi-exponential in time, reflecting two distinct eigenvalues) is reminiscent of Motulsky-Mahan competition dynamics at a monomer (Motulsky and Mahan 1984).
In this paper, we develop mathematical models for the dynamics of ligand binding at pre-dimerised receptors. In Sect. 2, we formulate and solve (analytically) a linear ODE model for single-ligand ( A) binding kinetics at constitutively dimerised receptors. We first relate our model to previous models and find an analytical equilibrium solution for total bound ligand. From this solution, we derive a condition under which multiple inflections in the logDR curve appear, in terms of the mechanistic binding cooperativity coefficient (α). Further, we show that the time-dependent problem has an analytical solution which may be easily constructed and computed without the need for numerical ODE solvers and use this solution to simulate the binding dynamics. In Sect. 3, we extend our model to account for the presence of a second competing ligand  B). Again, we relate to previous models and begin by finding the equilibrium solution. This enables us to derive a condition for multiple inflections in the logDR curve for total bound ligand A, which this time depends on multiple mechanistic binding cooperativity coefficients and the concentration of B. The ODE model for competition dynamics is linear as before, but its analytical solution is more laborious to compute. We again simulate the binding dynamics to explore the effects of the dynamic cooperativity factors. In Sect. 4, we conclude with a discussion of our main results, underlining our contribution to the literature, and considering dynamic modelling in the context of mathematics supporting drug discovery.

Model Formulation
In this initial model, we consider a single-ligand, A, binding to a pre-dimerised receptor complex. For the purposes of this paper, we assume all dimers are homodimers and so consist of two identical receptors bound together. We denote these dimers as R for simplicity in the model but keep the assumption that ligand molecules can bind to either side of the dimer. Association and dissociation rates are k a+ and k a− , respectively, and we assume there is no bias to which protomer the ligand molecule binds to first. As the second ligand binding may be affected by one side of the dimer being already bound, we introduce α = α + /α − which represents the equilibrium binding cooperativity, that is the factor change in affinity for the dimer when it is already ligand bound. The value α = 1 represents neutral cooperativity, and α > 1 and α < 1 represent positive and negative cooperativity, respectively. Figure 1 shows the system of biochemical reactions. Since R represents a dimer, AR is the complex created by a single-ligand molecule binding to a protomer and AR A is a dual bound receptor.

Differential Equations
Applying the law of mass action, the binding kinetics are governed by the following system of ordinary differential equations (ODEs): with initial conditions

Equilibrium Analysis
In the spirit of classical receptor theory, we first investigate the equilibrium behaviour of the system, in particular the effect of the equilibrium cooperativity factor α. The equilibrium relationships are where K A = k a+ /k a− is the equilibrium association constant and α = α + /α − is the equilibrium binding cooperativity. The total concentration of dimers is which we can combine with equations (4) to express the equilibrium concentrations in terms of parameters only, giving This result is equivalent to the result (Giraldo 2008) states when discussing a mechanistic model for a single-ligand binding to a dimerised receptor. As the maximal ligand bound is 2D tot , we can calculate the EC 50 value for A bound (the concentration giving half-maximal effect) as (8) Figure 2 shows how the equilibrium cooperativity factor α affects the log doseresponse (logDR) binding relationship. Increasing cooperativity leads to dimers becoming dual bound for lower concentrations of [A] and lower peaks in [AR]. Further, for A bound , biphasic logDR curves with multiple inflections are possible. With extreme negative cooperativity, we see three inflections in the curve instead of just one. It can be shown ("Appendix B" ) that the inflection point that appears in all curves is at the point [A] = 1/K A √ α, which in this case is the same value as the EC 50 due to the symmetric nature of the curves. Also, with increased α > 1, we see an approximate leftward shift in the curve.
Since the existence of extra inflections depends on α we seek a condition on α for when these appear. Following the calculation in "Appendix B," we find that for three inflections we require We note that this is different from Giraldo's dimer cooperativity index condition for apparent cooperativity (Giraldo 2008). From a curve-fitting perspective, α = 1/4 is mechanistically negatively cooperative, but would give apparent cooperativity in the Hill function sense (see Giraldo 2008). When (9) holds, the inflection at [A] = A 50 changes from a rising inflection point to a falling one and we get two extra inflection points, one for [A] < A 50 and one for [A] > A 50 at the points We see in Fig. 2 that this results in a biphasic curve. We note here that a biphasic logDR curve with three inflections as shown would not be seen for monomeric receptors. An experimental logDR curve may be suggestive of pre-dimerised receptors and very negative cooperativity.

Binding Dynamics-Analytical Solutions
The exact solution to the initial value problems 2, 3 may be constructed by calculating the eigenvalues of the coefficient matrix and by using the method of undetermined coefficients ("Appendix D"), or by the method of Laplace transform (as in Motulsky and Mahan 1984). Both methods give the analytical solutions as: where the eigenvalues of the coefficient matrix M are and Tr(M) and det(M) are the trace and determinant of M, respectively, given by: Our response of primary interest is the total concentration of ligand bound, A bound , which is given by The solution has two exponential components, meaning we would expect to see biexponential time courses in general. It can be shown ("Appendix C") that not only are eigenvalues λ 1 and λ 2 real and distinct but they are also always negative, since all parameters in the model are positive. Further, λ 1 gives the slow component of the dynamics, while λ 2 gives the fast component. Clearly, as t → ∞, the exponential components will decay to zero and we recover the equilibrium concentrations as giving the total ligand bound at equilibrium as equivalent to (6). Having analytical solutions for the dynamics allows time courses to be plotted without the need for numerical ODE solvers. This is particularly useful for pharmacologists without numerics expertise, allowing them to construct exact solutions for any parameter values in software packages such as Microsoft Excel or GraphPad Prism (accepted). In this respect, the dynamics for our dimer binding model have the same analytical structure as Motulsky-Mahan dynamics (Motulsky and Mahan 1984). In "Appendix H," we show a sample MATLAB code for construction of the analytical solution.

Results
Here we present numerical results which demonstrate the effects of cooperative dynamics across dimerised receptors. Figure 3 shows the binding of ligand A to the dimerised receptor [R] with the ligand being added as a constant quantity. Initially, ligand molecules bind singly to dimers, creating an increase in [AR]. As time increases, a positive cooperativity factor means that the chance of a second ligand molecule binding to the singularly bound dimer is increased; thus, [AR A] increases. This increase in [AR A] also has the effect of reducing [AR] which then falls towards zero. Hence, we see a peak in [AR]. In fact, this peak becomes a point of interest as we move on For a peak in [AR](t) (see "Appendix E"), we require and the corresponding peak concentration of [AR] is Binding dynamics with varying α. As we move from positive to negative cooperativity, we see a more pronounced peak in [AR] with A bound tending to lower concentrations. Here, In Fig. 5, we study the effect of α on A bound (t). Clearly, the binding rate for a second ligand to the dimer increases and decreases for positive and negative cooperativity, respectively. We see a peak in [AR](t) whose timing and concentration are dependent on α + . As α + increases, the peak in [AR] decreases and occurs earlier. Both [AR A] and A bound increase as α + increases.

Model Formulation
We now introduce a second ligand, B, to the system. The rationale for this is that quantification of effects of unlabelled ligands can be achieved by competition experiments with labelled and unlabelled ligands, as seen in previous studies involving monomeric and dimeric receptors (Motulsky and Mahan 1984;May et al. 2011;Durroux 2005;Franco et al. 2006). The kinetics of this system are key in highlighting and quantifying allosteric interactions across dimerised receptors, as indicated by May et al. (2011), who discuss the influence of an unlabelled ligand on the dissociation (washout) kinetics of a labelled ligand, when dimers are present.
We denote β = β + /β − as the influence a protomer bound by ligand B has on a second B molecule binding, and γ = γ + /γ − as the cooperativity factor describing the interaction between A and B bound receptors. This extended system gives a set of six reactions, as can be seen in Fig. 6. Taken together, these reactions correspond to a single face of the model schematic developed by Franco et al. (2006).
Schematic showing the binding possibilities with two ligands in the system

Differential Equations
As the total concentration of receptors is given by we may eliminate [R] to give a system of five linear ODEs: We note that this system is a simplification of the one developed by May et al. (2011) (see "Appendix A").

Equilibrium Analysis
Once the system reaches equilibrium, we have the set of relations The receptor conservation law therefore gives that Here we consider, as in May et al. (2011), competition experiments where ligand A is labelled and ligand B is unlabelled. We can clearly see that, providing [B] is fixed, The total concentration of bound ligand A (assumed an experimentally measurable quantity) is which, using Eqs. (24) and (25), gives As the maximal A bound for varying [A] remains as 2D tot , we calculate the EC 50 for A bound as which we note is independent of the A-B cooperativity factor γ . In Figs Fig. 10).
In Fig. 9, we show the effects of A-B cooperativity γ and see that, for extra inflections, γ is required to be high as opposed to low. We again require there to be a high concentration of B in the system. The individual species curves in Fig. 11  Investigating these extra inflections in A bound ("Appendix G"), we find that there is always an inflection point at We further find that extra inflections appear under the condition Again we have an extra inflection at either side of the original one, at the points

Time Course Results
Analytical solutions for the system (22) are theoretically possible, given that the system is linear. However, the task of computing eigenvalues exactly becomes laborious and impractical; we instead construct the solutions by numerical evaluation of the eigenvalues. A numerical ODE solver could also be used. Here we present time course simulations for two-ligand competition. In Fig. 12 In Fig. 13, we see the effect α has on the time course dynamics of the system. We look at a range of values of α while keeping β and γ fixed, giving neutral cooperativity.

It is clear that regardless of the A-A cooperativity, we get peaks in both [AR] and [B R].
As α increases, the peak in [B R] is lower and the peak in [AR] is later and higher. As In Fig. 14, we vary cooperativity factor β. Non-monotonic behaviour is clearly shown for all receptor states except [B R B], and the potential for non-monotonic signal A bound is clearly demonstrated.

Discussion
In this paper, we have presented dynamic models of ligand binding to pre-dimerised GPCR homodimers, for both a single-ligand and two-ligand competition. The models are linear ODE systems, allowing analytical solutions for time-dependent and equilibrium responses. The model formulation, solution and results serve as a contribution to the field of pharmacological modelling and are expected to be of practical use, given the ease with which we can compute solutions. In particular, the bi-exponential single-ligand binding kinetics show a similar model solution structure to the widely used Motulsky-Mahan model for competition binding at monomers (Motulsky and Mahan 1984). We therefore propose that this model can be adopted, interpreted and implemented with relative ease by pharmacologists, and as such, we have provided a recipe for computational solution.

so cooperativity is neutral for A-A and B-B and A-B cooperativity depends solely on γ +
Time course computation for our models is straightforward, given the analytical solutions we have developed. However, there are interesting features in the equilibrium logDR curves for the experimental readout (signal) A bound . We have noted the possibility for multiple inflections in logDR curves, for both single-and two-ligand binding. Multiphasic features in logDR curves are not reproducible by standard Hill functions (Veroli et al. 2015) which only support single inflection curves. For singleligand binding, multiphasic logDR curves theoretically rule out monomeric receptor binding as the only ligand binding mechanism; such experimental data therefore would suggest another binding setup, possibly due to dimeric receptors. The existence of multiple inflections depends on the level of cooperativity across a dimerised receptor, and importantly we are able to give a condition on the single-ligand cooperativity factor for a three-inflection curve (9). The practical use of this condition is clearly in assessing the sign and magnitude of cooperativity towards quantitative classification of drug-receptor interactions. Given that Hill functions are not suitable for fully characterising multiphasic logDR curves (Veroli et al. 2015), our analysis here goes beyond dimer cooperativity indices which stem in part from cooperativity in the Hill function sense (Giraldo 2008); the present work concerns mechanistic cooperativity which is explicit in the original model schematic, as opposed to empirical measures from a more limited model. In Veroli et al. (2015), it is noted that multiphasic logDR curves have particular importance in a number of contexts including cancer pharmacology, and effort should be made to move beyond Hill function fitting wherever possible.
Two-ligand competition at GPCR homodimers has previously been discussed and modelled (May et al. 2011). In addition to the therapeutic relevance of agonistantagonist competition (May et al. 2011;Bridge et al. 2010), for homodimeric receptors, dissociation kinetic assays may represent a more sensitive forum for detecting negative cooperativity than logDR analysis alone. As such, we have reduced the model of May et al. (2011) and outlined the analytical solution structure for the time-dependent problem. Further, considering a system with a labelled ligand A and an unlabelled ligand B, we are able to find conditions on [B] and the three cross-dimer cooperativity factors for multiple inflections in the logDR curve for A bound (Eq. 31), which again is designed to aid interpretation of experimental data and clear identification of parameter regimes which may have therapeutic significance.
Our simple models of ligands binding to constitutive homodimers, and particularly the explicit statement and elucidation of analytical results, represent an important addition to receptor theory. We remark that this should be considered only the foundation of the theory for binding and activation dynamics of GPCR dimers. There is a clear pathway for building on this foundation to expand upon this model towards simulating and understanding the potential diversity binding and signalling outcomes of non-monomeric receptors. In particular, the following natural extensions and considerations should be explored in future work.
1. Heterodimeric GPCRs are widely acknowledged to exist (Smith and Milligan 2010;Milligan 2004Milligan , 2006Milligan , 2009. Here, we have reduced the model of May et al. (2011) using the symmetry inherent in their homodimer model, resulting in a clean formulation amenable to algebraic manipulation. Extension to a more general model which will encompass heterodimers will require reversion to the May et al. schematic, but with an implicit asymmetry which will allow for an "extra level of complexity" (Bai 2004).
2. Higher-order complexes may exist (Smith and Milligan 2010;Juška 2008), and while a general functional signalling model for an nth-order oligomer may be impractical, a ligand binding model of the type developed in the current work should be possible and may be valuable in distinguishing between orders of oligomer. 3. Here, we have assumed that all dimeric receptors are constitutively produced, are the only receptor entities and have not considered the dimerisation process itself. However, a mixture of monomers, dimers and higher oligomers may exist (Smith and Milligan 2010;Milligan 2004), and ligand-induced effects on the dimerisation process have also been noted for GPCRs (Smith and Milligan 2010;Bai 2004) and other receptors (Mac Gabhann and Popel 2007). Therefore, we propose to expand the model to include constitutive and ligand-induced dimerisation dynamics. This will allow theoretical investigation of questions put forth in Milligan (2004) regarding the effects of fraction of receptors which are dimers, dynamics and lifetime of dimerisation and the effect of ligands. Further modelling in this direction will draw on dimerisation models (Vera et al. 2008;Mac Gabhann and Popel 2007), towards a general, hybrid model for constitutive and ligand-induced dimers and ligand binding. Further complexity may be possible within a dimer population, whereby only a fraction of dimers may support cross-dimer crosstalk (Durroux 2005). 4. In keeping with classical analytical pharmacology approaches, we have first considered models for ligand binding; the cooperativity factors α, β, γ are binding cooperativities. A diverse range of pharmacological effects is theoretically possible from dimeric receptors, in terms of receptor activation and signalling dynamics (Milligan 2004).
Much is yet to be understood about signalling downstream of dimerised receptors, and how to distinguish between ligand-receptor level crosstalk and downstream crosstalk when interpreting experimental data (Milligan 2006(Milligan , 2009Chabre et al. 2009). A natural extension of our model will include receptor activation (Giraldo 2008), G protein binding (Milligan 2009;Bai 2004), and ultimately G protein activation and cycling (Woodroffe et al. 2010;Bridge 2009). In moving towards simulating downstream functional effects, new modelling studies will require consideration of cross-dimer cooperativity in (i) ligand binding, (ii) receptor activation and (iii) G protein binding. 5. Beyond simulating and using our theoretical results to identify parameter ranges, our models will be useful for parameter estimation to return best-fit parameter values to experimental data, as in May et al. (2011). Considering time course data for ligand binding, and possibly washout experiments, it is important to determine which of the model's unknown kinetic parameters are theoretically identifiable from the given readout. In addition to model fitting with any available data and pseudo-experimental data, using methods such as genetic algorithms (Ashyraliyev et al. 2009), an important (but often overlooked) computation in bio-modelling is that of structural identifiability analysis (SIA). We propose to investigate the identifiability properties of our linear models using Laplace transform methods (Godfrey 1986). While beyond our scope here, this is an important avenue of investigation, given that in May et al. (2011), an initial parameterisation of the model therein suggests a non-identifiability issue; estimated equilibrium cooper-  (May et al. 2011) ativities are reported rather than all individual kinetic constants. The analysis in May et al. (2011) further suggests that a combination of association and washout experiments will be useful in estimating kinetic parameters relating to dimerised receptor binding. The simplified models presented in the current work will be used to simulate such experiments, and a detailed SIA will be performed towards parameterising ligand-dimer interactions. In Fig. 16, we show an initial simulation of binding followed by washout, in qualitative agreement with the results presented in May et al. (2011).
The current work provides a theoretical foundation for the study of ligands binding dimerised receptors. We consider this theory as a vital step towards the potential exploitation of dimerised receptors in a therapeutic and drug discovery context. Understanding the functional, physiological and therapeutic significance of receptor dimers is an ongoing challenge (Smith and Milligan 2010), with some questions still outstanding over the possible requirement of dimerisation for receptor function (Milligan 2004). Given the broad appreciation of GPCR dimerisation, it is now widely accepted that dimers are true drug receptors, and cross-dimer binding and activation cooperativity gives rise to a "new kind of pharmacological target" (Ferré et al. 2014) whose allosteric cooperativity is key. So-called dual molecules can be used in targeting dopamine-adenosine heterodimers in the treatment of Parkinson's disease, while it is suggested that opioid-cannabinoid heterodimers could be targeted for pain relief (Franco et al. 2008). Mathematical modelling will undoubtedly continue to be a key tool in simulating and quantifying in vitro ligand-receptor interactions and towards understanding of in vivo dynamics. In view of the potential functional diversity of dimerised receptors, the further development of dimerisation dynamics as a field of receptor theory is a crucial mathematical challenge towards drug discovery. Acknowledgements CW is supported by a Swansea University/College of Science PhD studentship. We acknowledge helpful discussions with Professor Venkat Kanamarlapudi and Dr. Vivi Rottschäfer.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Simplification of the May et al. Model
System (22) is seen to be a simplification of the model developed by May et al. (2011), equivalent by taking noting that where a tilde denotes a parameter or variable in the May et al. (2011) model.

B Extra Inflections in the Single-Ligand A bound Curve
In Fig. 2, we saw extra inflections in the logDR curve for A bound . We investigate these and derive the condition given in Eq. (9). Writing A bound in terms of X = K A [A] for simplification, we have To investigate the inflections, we require (and X = K A 10 Y ). Using Di Bruno's Johnson (2002) formula for the chain rule for second derivatives, we find Calculating each term individually, we have dX dY = K A 10 Y ln 10 = X ln 10, Substituting into (38), and after some simplification, we find To locate the possible inflection points, we find zeros of d 2 A bound /dY 2 , which occur when This gives the trivial solution at X = 0 and solutions at However, we are only interested in positive solutions, so we conclude that we always have an inflection point at X = 1/ √ α, that is when We note that this is the same concentration of [A] that gives half the maximal response. Other possible inflections occur for roots X > 0 of the quadratic equation which gives solutions at We see that we get a repeated root if that is when α = 1/16 or α = 1/4. However, we see that the root is only positive for α = 1/16. To get real and distinct solutions, we require that is when α > 1 4 or when α < 1 16 . We also need to confirm whether the roots are positive or negative under these conditions. So first noting that Now if α < 1/16, then 1 − 8α > 1 − 1/2 = 1/2 > 0; thus, we have two positive roots under this condition. However, if α > 1/4, then 1 − 8α < 1 − 2 = − 1 < 0; hence, under this condition we get two negative roots. So we conclude that we get two possible extra inflection points under the condition 0 < α < 1/16. To confirm that these are inflection points and classify the type of inflection, we consider d 3 A bound /dY 3 . Using the Di Bruno formula (see Johnson (2002)) for higher derivatives, or MATLAB's Symbolic Toolbox, we find: For the EC 50 value, at X = 1/ √ α, and substituting this into the third derivative. After some simplification, we have which we can see is positive for any 0 < α < 1/16 and negative for any α > 1/16. Similarly for the extra inflections, we found at X = τ 1 and X = τ 2 when α < 1/16. Substituting both of these in results in a value of (ln 10) 3 (16α − 1) which clearly gives a negative result for any 0 < α < 1/16. Looking at these results, we can see that if α > 1/16, we have the single inflection point, which in this case has a concavity that changes from convex to concave. However, when 0 < α < 1/16, this inflection point reverses and so is concave changing to convex. In this case, we also get the two extra points, with each of those being convex changing to concave at the point.

C Confirming the Eigenvalues of the Coefficient Matrix are Real, Distinct and Negative
Before we can solve the system of equations to give analytical solutions, we need to know whether the eigenvalues, given in Eq. (13), of the coefficient matrix M are real and distinct so we can determine the form of the general solution. Being able to state whether these are also positive or negative allows us to analyse these solutions. The eigenvalues are defined as so we can see that the discriminant is Expanding D, we have Hence, D is always positive, providing that all parameters are positive.
We can also show that both λ 1 and λ 2 are always negative. First noting that Tr(M) is negative due to all terms being negative, and

D Analytical Solutions for the Single-Ligand Dynamics
One of the main motivations for this paper is to be able to give exact solutions to the system, allowing for simulations to be run easily; hence, in this appendix we derive these solutions as stated in Eqs. (11) and (12). To find a solution to the system, we first find the eigenvalues of the coefficient matrix, M.
The eigenvectors corresponding to eigenvalues λ 1 and λ 2 are Thus, we can state that a solution to the system will be of the form where X p (t) is a particular solution, found by solving MX + f = 0, thus As t → ∞, all other terms in the solution tend to zero, due to λ 1,2 always being negative; thus, the particular solution tells us what [AR] and [AR A] will be at equilibrium. We now use the initial conditions to solve the system and find the coefficients c 1 and c 2 , as We find the solutions

E A Condition for Existence of a Peak in [ AR](t)
As we saw in Figs. 3, 4, we sometimes get a peak in [AR](t). In this appendix, we investigate these peaks to find a condition Eq. (19) for peak existence and calculate this peak concentration Eq. (20) is. The peak time is given by Eq. (18) For peak existence, we therefore require where λ 2 < λ 1 < 0. Hence, we require, for a positive peak time, which requires that α − k a− + λ 1 < 0. Since , where D is the discriminant for the coefficient matrix, we find our requirement to be which eventually gives To evaluate the peak concentration [AR], we take Eq. (18) and substituting the peak time into Eq. (11), giving as the corresponding concentration of [AR].

F Individual Species Dose-Response Curves
To further help us understand the behaviour observed in the logDR curves for A bound , we look to the individual species plots. In Fig. 17, we see plateaus in the logDR curves of [AR] and [AR A] when we have very low A-A cooperativity. These are the main contributions to the inflections in A bound .

G Investigating the Extra Inflections in the Two-Ligand A bound curve
Similarly to the single-ligand system, we also get extra inflections in the A bound curve when we have two ligands in the system. In this appendix, we investigate these to derive the condition given in Eq. (31) under which we get these extra inflections. We begin by writing A bound as = X (ln 10) 2 (n − α X 2 )(mn + (8αn − m 2 )X + αm X 2 ) (n + m X + α X 2 ) 3 where m = 1 + γ Y and n = 1 + Y + βY 2 . Clearly, we have the trivial root at X = 0 as well as a root when n − α X 2 = 0, that is when As we are only interested in positive roots, we conclude that we have a possible inflection at which corresponds to the EC 50 concentration. We also get extra roots when mn + (8αn − m 2 )X + αm X 2 = 0, which using the quadratic formula gives us X ± = m 2 − 8αn ± (8α − m 2 ) 2 − 4αm 2 n 2αm (85) = m 2 − 8αn ± (m 2 − 16αn)(m 2 − 4αn) 2αm .
We see that we get a single, real root when that is when m 2 = 16αn or m 2 = 4αn. However, we see that if m 2 = 4αn, then X = − 2n/m, hence a negative root. Whereas if m 2 = 16αn, then X = 4n/m, so we have a positive root. Thus, we can conclude that we have a possible inflection at X = 4n/m under the condition of m 2 = 16αn. Finally, we have two real distinct roots when (m 2 − 16αn)(m 2 − 4αn) > 0.
Using this, we can state that if m 2 > 16αn, then m 2 −8αn > 16αn −8αn = 8αn > 0; hence, we have two positive roots. Conversely, if m 2 < 4αn, then m 2 − 8αn < 4αn − 8αn = − 4αn < 0; hence, both roots are negative. Thus, we can conclude that we get two positive roots, and possible inflection points, under the condition m 2 > 16αn. That is To confirm that these are in fact inflection points, we first confirm that X − < X c < X + . First noting that m 2 −16αn < m 2 −4αn ⇒ (m 2 −16αn) 2 < (m 2 −16αn)(m 2 − 4αn) ⇒ m 2 − 16αn < (m 2 − 16αn)(m 2 − 4αn), then we can state that = 16n 2 m 2 (97) Thus, we can confirm that X − < X c . Also hence confirming that X c < X + . We now use the second derivative to confirm there is a sign change at each of the points. We can see as X → 0 + , we have d 2 A bound /d Z 2 > 0. TakingX − = 4n/m as a point between X − and X c , we evaluate d 2 A bound dZ 2 X =X − = − 12n 3 (m 2 − 16αn) 2 m 4 < 0; thus, we can confirm there is an inflection point at X − with the curve changing from convex to concave. TakingX + = (m 2 − 8αn)/2αm as a point between X c and X + , we evaluate d 2 A bound dZ 2 X =X − = (m 2 − 8αn)((m 2 − 8αn) 2 − 4αm 2 n) 2 32α 3 m 4 > 0; hence, there is a sign change at the point X c with the curve changing from convex to concave. Finally, taking the limit as X → ∞, we see that d 2 A bound /d Z 2 < 0; thus, we have a third inflection point with the curve changing from convex to concave. So to conclude, under the conditions 16α(1 + Y + βY 2 ) < (1 + γ Y ) 2 , we have three inflection points.