Classical structural identifiability methodology applied to low-dimensional dynamic systems in receptor theory

Mathematical modelling has become a key tool in pharmacological analysis, towards understanding dynamics of cell signalling and quantifying ligand-receptor interactions. Ordinary differential equation (ODE) models in receptor theory may be used to parameterise such interactions using timecourse data, but attention needs to be paid to the theoretical identifiability of the parameters of interest. Identifiability analysis is an often overlooked step in many bio-modelling works. In this paper we introduce structural identifiability analysis (SIA) to the field of receptor theory by applying three classical SIA methods (transfer function, Taylor Series and similarity transformation) to ligand-receptor binding models of biological importance (single ligand and Motulsky-Mahan competition binding at monomers, and a recently presented model of a single ligand binding at receptor dimers). New results are obtained which indicate the identifiable parameters for a single timecourse for Motulsky-Mahan binding and dimerised receptor binding. Importantly, we further consider combinations of experiments which may be performed to overcome issues of non-identifiability, to ensure the practical applicability of the work. The three SIA methods are demonstrated through a tutorial-style approach, using detailed calculations, which show the methods to be tractable for the low-dimensional ODE models. Supplementary Information The online version contains supplementary material available at 10.1007/s10928-023-09870-y.


Introduction
Mathematical models of pharmacological systems have become key in understanding the interactions between ligands and living cells, and as such play a significant role in the development of new therapeutic medicines.These models are often comprised of ordinary differential equations (ODEs) which depend on mechanistic parameters that represent biological processes and whose values are largely unknown [13].An essential step in using these models requires estimating values of these parameters [52] by fitting to experimental data from, for example, ligand binding assays.Parameter estimation for ODE models of biological systems typically involves optimisation algorithms [17].However, these fitting routines can result in inaccurate and unreliable estimates [17,42].
Identifiability analysis is the process of assessing whether it is even theoretically possible to estimate a set of parameters uniquely from experimental observations and the dynamic equations [20,52].Such analysis is therefore required to determine the reliability of parameter estimates.In particular, structural identifiability analysis (SIA) uses the model structure, together with observed outputs, to determine whether parameters can be returned successfully, given perfect, noiseless and bias-free observations [3,51].

A simple demonstration of the issue of (non-) identifiability
Theory and methodology for identifiability analysis have been developed and considered within the context of compartmental models for pharmacokinetics (PK) [4,14,22,23,29,30].We use a simple PK model to demonstrate the issue of non-identifiability.Consider the schematic shown in Fig. 1, which depicts a two-compartment model for oral absorption of a drug.In response to a drug input rate u(t) (eg. a bolus dose), the drug amount in the absorption compartment is a b ðtÞ and the drug amount in the central compartment is a c ðtÞ.
Assuming linear pharmacokinetics gives the following linear system of ordinary differential equations (ODEs) for a b ðtÞ and a c ðtÞ: where u(t) is the input dosing rate.It is straightforward to show [22,27] that the drug level in the absorption compartment following a single bolus dose D 0 (so that uðtÞ ¼ FD 0 dðtÞ, where F is the bioavailability factor [22,27,31] and d is the Dirac delta function) is given by a b ðtÞ ¼ D 0 Fe Àk a t : Note that an equivalent problem is given by setting uðtÞ 0 and setting a b ð0Þ ¼ FD 0 .The drug level in the central compartment is found to be a c ðtÞ ¼ D 0 Fk a k a À k e e Àk e t À e Àk a t À Á : In PK studies using such models, a c ðtÞ (or some scaled drug amount, for example, drug concentration given by a c =V where V is the volume of distribution [27]) would be of primary interest and would be the measured/observed output of the ODE system.Suppose a dose D 0 ¼ 500 (in arbitrary units) is given at time t ¼ 0. In Fig. 2, we show a c ðtÞ and a b ðtÞ for the two different sets of parameter values (k a , k e and F) given in Table 1.Despite the different parameter values, the two parameter sets give identical observed output a c ðtÞ (despite the unobserved, intermediate timecourses for a b ðtÞ being different).Hence the three parameters k a , k e and F are theoretically unidentifiable from such a time course; their values cannot be determined without further a priori knowledge.This system in fact gives an example of local identifiability, where multiple distinct values, unique in a neighbourhood, solve the problem -see definitions, e.g., in [12].Clearly, issues of identifiability should be considered as part of any parameter estimation implementation using real data.

Classical structural identifiability analysis (SIA) methods
The origins of SIA lie in the works of Kalman [33] for linear ODE systems, and Hermann and Krener [25] for nonlinear models.Since these works, many methods have emerged to assess the identifiability of a given model.For linear models, the Laplace transform (transfer function) method [4] may be used.For linear or nonlinear models, Taylor series expansions [45] and similarity transforms [8,28] are applicable.In theory, the methodology underlying these three approaches applies to an n-dimensional ODE system.In practice, the complexity of the algebraic manipulation involved in the required computations may limit their applicability to low-dimensional systems.Often compartment models arising in pharmacokinetics are used as suitable low-dimensional examples of the application of SIA methods [10,21,50].SIA methodology has grown as a field, and a number of methods and software packages have been developed, which are particularly useful for nonlinear or high-dimensional systems, many making use of symbolic algebra [2,5,12,13,[37][38][39].More recent SIA algorithms include singular value decomposition of sensitivity matrices [32,48] and scaling approaches [7].While the field of SIA has grown, the difficulties in applying both classical and newly developed methods in general are still noted throughout the literature.The intractability of the required algebraic manipulations for large and/or nonlinear ODE systems has proven to be a ''persistent bottle-neck'' [48] and most of the available methods may be ''too complex mathematically for the general practitioner'' [7].It is this limitation of SIA methods that has resulted in only a very small proportion of theoretical biology studies considering the identifiability issue [7].
Receptor theory is a core component of pharmacological analysis which considers the interactions between ligands and receptors, and the implications of these interactions, at the top of signal transduction pathways of pharmacological interest.A key aim of analytical pharmacology is to use timecourse data to estimate kinetic association and dissociation parameters to quantify ligand-receptor interactions.Similar to in many systems biology modelling work, SIA is frequently overlooked for ligand binding models.The aim of the current work is an informative tutorial which brings classical SIA methods to the field of receptor theory, in particular considering the dynamics of ligand binding models of importance.

u(t) a b a c k a k e
Fig. 1 A two-compartment PK model for oral absorption of a drug SIA has been applied to pharmacodynamics models with drug effect as the observed output [31], and to more complex PK/effect models with measured outputs beyond/downstream of ligand binding [11,16].Ligand binding has been considered in a PK context with total ligand (free ligand plus bound ligand) as the output [9].Identifiability for equilibrium ligand binding is only considered in detail in [42,43].Here we follow a receptor theory approach and consider SIA for models of ligand binding assay dynamics with total ligand-bound receptor as the readout for binding scenarios of practical importance [34,44,55].Given the similar size and structure of many ligand binding receptor theory ODE models [34,36,54] to the PK models for which the three ''classical'' SIA methods have been demonstrated, we propose that application of these classical methods to widely used receptor theory models is valuable on two fronts.Firstly, the analysis will provide new results on the identifiability of key kinetic parameters.Secondly, the study will bring the SIA methodology to a new audience in receptor theory.

Paper overview
Here we consider three ligand-binding scenarios.The first concerns a single ligand type binding monomeric receptors, which is the starting point for much of receptor theory [34,35].The second scenario is the widely used competition binding model of Motulsky-Mahan [44], wherein two different ligand types compete for monomeric receptors.This model can be formulated as a second order linear ODE system with four kinetic rate constants and the total receptor concentration as the ''unknown'' parameters.It is known intuitively that a subset of the parameters may be estimated if the other parameters are already known [44], but formal analysis of what is theoretically identifiable from a single experiment has not been presented.Furthermore, due to its widespread use, the practicality of parameter estimation for this model continues to receive attention [15,18,49].Theoretical (SIA) questions still remain.
The third scenario we consider is that of a single type of ligand binding to homo-dimerised G protein-coupled receptors (GPCRs).This model is a natural extension of the monomeric receptor model to account for the growing acceptance of the existence of GPCR dimers.This model may also be formulated as a second order linear ODE system.An analysis of the binding kinetics appears in [55], but SIA questions have not yet been addressed.
The remainder of this paper is organised as follows.In Sect.2, we introduce the three ''classical'' methods by their application to the simple model for ligand binding to monomeric receptor.In Sect.3, we apply the three classical methods to the Motulsky-Mahan and GPCR dimer models, yielding new identifiability results of practical importance to the pharma-modelling community, and also contrasting the three approaches.In Sect.4, we extend the analysis to consider further experiments (including equilibrium binding, washout and multiple timecourses) which can be performed to overcome issues of non-identifiability.We conclude in Sect. 5 with a discussion of our main results and contribution to the literature.

Methods: applying SIA to monomeric receptor binding model
To demonstrate the methods, we use a simple model of a monomeric receptor binding with a ligand.While the analysis of this model is mostly trivial, it will allow us to    explain the process of using each method to determine identifiability.Furthermore, this model has been the basis of much of the theoretical foundations of receptor theory of the past few decades and is still widely used in drug development research [34,35].We assume that a ligand, A, binds to the monomeric receptor R with association and dissociation rates k aþ and k aÀ respectively, to create the complex AR.This can be described as a chemical reaction as:

AR
The law of mass action, assuming ½A is held constant, gives the linear system of ODEs Now, asssuming that the only known value is These ODEs, together with the initial conditions form the initial value problem describing the kinetics of the system, where R tot is the total receptor concentration.Throughout our work, we focus on constant ligand concentration, linear models, in line with typical, classical analyses in the literature ( [34,44]).We note that the lesscommon assumption of of significant ligand depletion effects would yield a nonlinear system for which the transfer function method (see below) is not applicable.
For the model 3, the concentration of ½AR is measured experimentally at a number of timepoints, hence the observed output is Our assumption here uses total bound ligand as a direct readout (e.g.[44,55]).More detailed models for specific binding assays could consider a further constant of proportionality or Hill function measured response [6,34].Now, asssuming that the only known value is the ligand concentration ½A, the vector of unknown parameters in the initial value problem (3a-d) is given by p ¼ ðk aþ ; k aÀ ; R tot Þ.
For this simple model, it is possible to solve the initial value problem and obtain exact solutions.The measured output solution is found to be From this we can see that there are two identifiable parameter combinations, namely Noting that the second of these expressions appears in the denominator of the first, we see that the identifiable parameter combinations may be simplified and listed as We would therefore expect identifiability methods to also find these two identifiable parameter combinations.In particular, we note that a single experiment does not show the primary parameters of interest, k aþ and k aÀ , as identifiable.

The transfer function method
The first method we apply is the transfer function method, which makes use of the Laplace transform.This method was first proposed by Bellman and A ˚stro ¨m [4] and is one that is simple in nature, but restricted to linear time-invariant systems.Before applying this method to the monomer-ligand binding model, we first note that, the only input to the system is captured by the initial conditions (3c).We choose to reformulate in such a way that the ligand input is captured by a forcing term in the ODE, as in similar analysis in PK and control theory [1].We use a conservation law to reduce the system dimensions, in line with our previous work [55].In (3a-d), total receptor is conserved such that which we use to reduce the system by eliminating ½R.
Upon substitution, we find the single differential equation with initial condition and output The system is now in the format The transfer function describing the input-output relation for the system P in (9) (see [1]), is where I is the n Â n identity matrix if x is an n-vector.For the model given by (10), we find that The parameter-dependent coefficients of powers of s (including the s 0 terms) in both the numerator and denominator of Q(s), give the uniquely identifiable parameter combinations.We find a vector of parameter combinations that are identifiable to be To determine identifiability we assume the existence of an alternate vector of parameters which would yield the same output as the original parameters p ¼ ðk aþ ; k aÀ ; R tot Þ. Denote this alternative vector by e p ¼ ð g k aþ ; g k aÀ ; g R tot Þ.
Then a second vector of parameter combinations is fðe pÞ.Now we set fðpÞ ¼ fðe pÞ and solve for p.Each parameter, p i , is deemed structurally globally identifiable if, in solving this system of equations, it is found that p i ¼ e p i .Similarly, the parameter is structurally locally identifiable if there is a fixed number of possibilities for p i , and unidentifiable if there are infinitely many possible solutions for p i .Clearly, in this example, we have three unknown parameters and only two identifiable combinations, and therefore, it is not possible to identify all parameters from a single output time course.In fact, solving fðpÞ ¼ fðe pÞ we find that none of the parameters are identifiable individually, and so we conclude that only the grouped parameters are identifiable.Non-identifiability for (8a-c) can be seen in the numerical results shown in Fig. 3 that were generated using the parameter values in Table 2.Here we see how, for example, three different parameter sets can result in the same measured output curve (½AR) but have different magnitudes of the free receptor concentrations (½R).In Table 2 we see that the values of the identifiable grouped parameters given in (14) are equal for all three parameter sets.In Sect.4, we discuss approaches we can take to compensate for the non-identifiability of individual parameters.

Taylor series method
The next method we present makes use of the Taylor series.The Taylor series method was first developed by Pohjanpalo [45], and can be applied to either linear or nonlinear systems.While the receptor theory models we present in this paper are linear, we will later analyse the identifiability properties of a nonlinear system for ligand binding [57].We note however, that while the Taylor series may be used to determine identifiability of nonlinear systems, the algebra involved in applying the method to these problems can be difficult [10].
To apply this method to our example model, we can use the system in the form given in (3a-d), or the reduced form (as in (8a-c)), with both giving the same results.We choose to use the full system (3a-d).The Taylor series approach exploits the fact that there is a unique Taylor expansion for a given output yðtÞ about t ¼ t 0 , and so the Taylor coefficients (we refer to these simply as coefficients throughout) give identifiable parameter combinations [13].It can also be shown that, for linear problems, the maximum number of coefficients needed to determine identifiability is defined as where n is the number of variables [13].Coefficients beyond the first ð2n À 1Þ terms in the Taylor series give no further information about identifiability.
In our example of monomer-ligand binding (3a-d) we have n ¼ 2, and will therefore need to calculate a maximum of three coefficients.Here, our initial conditions are taken at time t ¼ 0, hence we take t 0 ¼ 0. Recall, the Taylor series of y about t ¼ 0, noting that we use the bracketed superscript to indicate the order of the derivative, is The first coefficient is simply As y ¼ ½AR, we use the initial conditions to obtain the first coefficient as Since the expression for yð0Þ contains no parameters, it gives no information about parameter identifiability.Since we have two remaining coefficients to determine but three unknown parameters, we can immediately conclude that the system is not globally identifiable.However, we will still continue to determine which, if any, parameters are identifiable and also any identifiable parameter combinations.The next coefficient, y ð1Þ ð0Þ, involves calculating the first derivative of the output function.Now, y ð1Þ ðtÞ ¼ ½AR ð1Þ ðtÞ: ð19Þ From (3b), we have Evaluating at t ¼ 0, using the initial conditions, then gives the second coefficient as Since ½A is known, we have where we denote c i as the i th identifiable parameter combination.Further derivatives are found with recursive substitution of the equations in system (3a-d).This gives the second derivative as After evaluation at t ¼ 0 and substitution of the initial conditions (3c), the second coefficient is found to be Using (22), we see that and since ½A is known, we may write as the second identifiable parameter combination.Hence, we have two identifiable parameter combinations.Moreover, these agree with the parameter combinations found when using the transfer function method, which is as we expected.To determine which parameters (if any) are identifiable, we would proceed as in the transfer function method by creating the vector fðpÞ ¼ ðc 1 ; c 2 Þ and solving fðpÞ ¼ fðe pÞ.
Calculating further derivatives gives no further identifiable parameter combinations.For example, the third derivative is, after some simplification, given by This gives the next Taylor coefficient as Fig. 3 Three sets of parameters are used to plot the system given in Eqs.(3a-d).All three parameter sets give the same measured output curve.However, non-identifiability can be seen in the individual species curves, ½R.Each set of plots is created using the values in Table 2 together with ½A ¼ 10 À8 M Table 2 The parameters for the three different parameter sets that are used to plot Fig. 3 Units Set 1 Set 2 Set 3 The identifiable parameter combination expressions in (14) are equal in each case and so is comprised of already identifiable parameters and parameter combinations.

Similarity transformation method
The similarity transformation method (or exhaustive modelling approach) was first proposed by Walter and Lecourtier [53], and originally was only applicable to linear problems.This was later extended to include nonlinear problems [50].The theory behind this method, along with proofs, can be found in [50].
In order to use the similarity transformation method, the system must be both controllable and observable [33].A system is said to be controllable if changing the input to the system changes the system states x, and observable if the initial state x 0 can be uniquely determined from a set of input-output measurements.In [50], a test for linear systems to determine whether the system is controllable and observable is outlined.For a linear system of the form ( 9), the controllability matrix C is defined by The system is then said to be controllable if rankðCÞ ¼ n, and observable if rankðOÞ ¼ n.
For our ligand binding model, we consider the system given by (9, 10), where the system input is captured in the ODE.That is, we have and n ¼ 1.As n ¼ 1, the system is both controllable and observable, hence we now move forward with checking identifiability.
The similarity transform method involves taking this system P ðpÞ, and assuming the existence of a system P ðe pÞ that depends on an alternate parameter set e p.These systems must be equivalent (that is, they give the same solution), and so must satisfy some equivalence conditions.The underlying theory and full algebraic equivalence theory can be found in [20,47] here we just state the conditions.We assume that there exists an n Â n matrix, T, that describes the transformation between the two systems.The conditions we impose are then det T 6 ¼ 0; ð32aÞ where a tilde indicates the alternate system.Applying these constraints on the two systems determines the identifiability of parameters.
In our system, we have n ¼ 1; and so T is a single entry matrix.As Clearly, det T 6 ¼ 0 and, as x 0 is independent of parameters, T e x 0 ¼ x 0 , and so conditions (32a) and (32b) both hold.
Applying condition (32d), we have Finally, we have, from condition (32c) Hence, we find the same identifiable parameter combinations as in the previous methods, with no individually identifiable parameters, confirming that all methods give the same results.

Results: SIA for further ligand binding models
Now that we have outlined the methods, we will apply each of these methods to the classical Motulsky-Mahan [44] competition binding model, in Sect.3.1, and to our recently presented dimeric receptor binding model in Sect.3.2.This will provide a tutorial of how each of the methods work in practice and also allows us to compare the different techniques.

Competition binding model
The Motulsky-Mahan [44] model is for a competitive binding scenario where two ligands, A and B, are each able to bind to a monomeric receptor.The ligand concentrations are assumed constant, and each ligand is able to bind the same receptor.Ligand A is radiolabelled and the total concentration of A bound to receptor can be measured experimentally.Ligand B is an unlabelled competitor whose total receptor-bound concentration cannot be measured directly.The reactions describing the interactions are: The system of ODEs governing the dynamics of the system is together with the initial conditions where R tot is the total receptor concentration.This forms the initial value problem describing the kinetics of the system.We assume that only the concentration of ½AR is measured experimentally, hence the output is The ligand concentrations ½A and ½B are known constants, and so a vector of unknown parameters in the model is p ¼ ðk aþ ; k aÀ ; k bþ ; k bÀ ; R tot Þ.In Sect.4.2, we will consider the scenario suggested in [44] where k aþ and k aÀ are already known from other experiments, but here we analyse the competition binding model with respect to the identifiability of all five parameters using a single timecourse.

Transfer function method
We first apply the transfer function method to determine identifiability of the parameters.In keeping with the Motulsky-Mahan analysis [44], we reduce the system in (36a-e) using the conservation law To write in the form (9), i.e., we introduce the matrices To check the identifiability of the system, we calculate the transfer function as which gives the vector of coefficients as We have five unknown parameters and only four identifiable combinations, therefore we can immediately conclude that it is not possible for all five parameters to be globally identifiable.Setting fðpÞ ¼ fðe pÞ allows us to determine identifiability of any individual parameters.The first two equations are Dividing ( 42) by (41) gives Hence, we find k bÀ to be identifiable, however, this is the only parameter to be so.This leaves the remaining identifiable parameter combinations as In Fig. 4 the non-identifiability of the model is clear, as three parameter sets (as given in Table 3) result in different species timecourse curves, yet the measured output curve of ½AR is the same for all sets.In particular, we notice that ½BRðtÞ has significant peaks in some of the curves, depending on the parameters used.

Taylor series method
To apply the Taylor series we consider the system (36a-e).
As we now have three state variables, we have n ¼ 3, so (15) gives that there is a maximum of five coefficients required to determine identifiability.Calculating the Taylor series coefficients, as stated in (16), and evaluating at t ¼ 0, gives the identifiable parameter combinations.Again, the first coefficient is which gives no information regarding identifiability.As we have five unknown parameters and only four remaining coefficients to evaluate, we can conclude that the system is not globally identifiable.The first derivative, as given by (36b), is giving the first Taylor coefficient as and as such, the first identifiable parameter combination as since ½A is known.Using recursive substitution of equations (36a-e), we can write the second derivative as Fig. 4 Three sets of parameters are used to plot the solutions of the system (36a-e).All three parameter sets give the same measured output curve, AR.However, non-identifiability can be seen in the individual species curves.Each set of plots is created using the values in Table 3 together with ½A ¼ 10 À8 M and ½B ¼ 10 À7 M 0.12 0.12 0.12 The parameter k bÀ , as well as the parameter combinations in the final three rows are equal in each case, as these are the identifiable quantities which gives the coefficient Now, and therefore, we have a second identifiable combination as Further coefficients are found in the same way, using recursive substitution of the system equations in (36a-e) to calculate higher order derivatives followed by substitution of the initial conditions.Using this method we obtain the third coefficient as which gives the third identifiable combination as We also obtain which, we find via some trial and error, may be written as which gives the identifiable parameter We note that the calculations are all performed using MATLAB Symbolic Toolbox [40].To conclude, we find that k bÀ is identifiable, and also the identifiable combinations Although these parameter combinations are not identical to those found from the transfer function method in expression (44), we note that hence the same four parameter combinations are indeed identifiable.

Similarity transformation method for competition binding model
In Appendix 1, we apply the similarity transform method, as introduced in Sect.2.3, to the Motulsky-Mahan competition binding system.We show that the same parameter combinations are found to be identifiable as for the previous methods (see (44) and ( 59)).
Comparing the methods, it is clear that, although all methods give the same identifiable parameters and parameter combinations, the transfer function method is by far the simplest in terms of ease of use.The Taylor series method, in particular, results in expressions that require quite some manipulation in order to obtain reduced expressions.

Pre-dimerised G protein-coupled receptor binding
The next model, and the final linear model, we consider is the GPCR homodimer model we presented and analysed in [55], for a single ligand binding.The schematic for the model is as follows.
Here, R represents the dimerised receptor, AR is the dimerised receptor with one ligand bound, and ARA is the dimerised receptor with both protomers bound by ligand.The parameters a þ and a À are the forwards and backwards binding cooperativities respectively.These capture the increased or decreased propensity for binding and dissociation when the opposite side of the dimer is ligand-bound rather than unoccupied.
The ODE system describing the model dynamics is given by (see [55]) with initial conditions The measured quantity is bound ligand, hence the output is given by We assume the only known parameter is the ligand concentration, ½A, and so we have the vector of unknown parameters as p ¼ ða þ ; a À ; k aþ ; k aÀ ; R tot Þ, where R tot is total dimerised receptor.In contrast to the monomeric receptor output (3d), we note that the output function is now a combination of two states, adding a significant difference to the proceeding computations.

Transfer function method
Again, we consider the transfer function method to determine identifiability.We first use conservation of receptors, which is given in this case by to reduce the system, giving This is in the form (9), where we identify We find the transfer function of the system to be Qðs; pÞ ¼ ð65Þ which gives the vector of identifiable parameter combinations as Hence we have no identifiable parameters but do have four identifiable parameter combinations.Again this can be seen in Fig. 5 where we show how three sets of different parameter values result in different individual species curves, yet all give the same measured output curve A bound .
While the curves of ½R and ½ARA are similar in shape across the three parameter sets, they have different magnitudes of concentration.The clearest differences are seen in the ½AR curves, where the different parameter sets result in curves that have distinctly different evolution patterns, with some curves having a peak and fall while others are monotonic.This highlights how naive parameter estimation performed without knowledge of identifiability issues could lead to incorrect conclusions being drawn about the underlying qualitative dynamics.The parameter values used for the plots are given in Table 4.

Taylor series method
We proceed with the Taylor series method to determine identifiability, with repeated substitution of the ODEs and initial conditions.While the process is the same as for the competition binding model in Sect.3.1, the output function being a combination of two state variables adds an extra complexity to the calculations.As in all previous sections, the first coefficient is trivial, that is The first derivative of the output function in (61e) is given by which, using the initial conditions in (61d), gives the first unique coefficient as The remaining three coefficients are calculated by following the method of Sect.3.1.2,by repeatedly differentiating the output expression and substituting in the dynamic equations (61a-e).We find that The parameter combinations in the final four rows are equal in each case; these are the identifiable combinations found in (66) A vector of identifiable combinations is then given by where While it is not immediately apparent that the Taylor Series method gives the same identifiable combinations as the transfer function method, we show in Appendix 2 that we can recover the four combinations in (66) from those in (73a-d) by algebraic manipulation, aided by symbolic computation.The result, again, is that the identifiable combinations are Similarity transformation method for GPCR dimer model In Appendix 3, we apply the similarity transformation method to the GPCR dimer model, and find that he identifiable parameter combinations are the same as those in (66) and (74).
Comparing the three methods applied to this system, we find that the transfer function method is the most straightforward to implement, whereas the Taylor series method results in expressions that require much simplification.

Results: addressing identifiability issues with equilibrium, washout and multiple time courses
The results thus far have shown none of the models to be globally identifiable from a single set of time course data.In this section we consider alternative ways in which all parameters can be identified.Commonly performed experiments for ligand binding include equilibrium (or saturation) binding assays, in which equilibrium binding levels are measured for a range of ligand concentrations to produce a concentration-response curve.For each ligand concentration, the binding experiments are run until equilibrium is assumed after which the amount of ligand bound is observed.These experiments are often used to estimate equilibrium constants K D ¼ 1=K A (the equilibrium dissociation constant), where K A ¼ k aþ =k aÀ , and R tot (total receptor concentration), for monomeric receptors [34].
Here, for each model in Sects.2, 3, we aim to establish identifiability for the corresponding equilibrium model, then use ''known'' equilibrium parameters together with timecourse data to establish identifiability for those kinetic parameters which were previously unidentifiable.
Washout experiments can also be used to gain further insights into the dissociation kinetics of ligands.In these experiments the free ligand is removed by repeated washing, ensuring that no further ligand associates with the receptors [41].Such experiments isolate the effect of dissociation and preclude further binding.Here, we consider this type of experimental dataset to establish identifiability of the kinetic dissociation parameters, in conjunction with association (binding) time course data to also determine identifiability of association parameters.
Finally, we also consider multiple binding experiments, whereby each data set is collected from experiments performed using different ligand concentrations.These data sets are then used simultaneously, with the aim of determining the minimum number of data sets required to make the model globally identifiable.In each case, we choose one identifiability method to apply.

Monomeric receptor binding with a single ligand
Recall the model for a monomeric receptor binding with a single ligand, as given by the schematic A þ R k aþ k aÀ AR with the system of equations as given in (3a-d).In Sect.2, we found that there are no identifiable parameters, only the parameter combinations are identifiable.We will use three different approaches to establish global structural identifiability, namely, concentration-response/saturation data together with a single set of time course data, a combination of association and dissociation data, and also multiple time courses.In each case we use an appropriate method from the three that we outlined in Sect. 2. In most cases this is the transfer function method due to its simplicity of implementation, however, dissociation timecourse data are analysed using the Taylor series method.

Equilibrium saturation curves
We first establish the identifiability of equilibrium parameters associated with ligand binding, namely K A ¼ 1=K D (the equilibrium dissociation constant) and R tot (total receptor).Note that, at equilibrium, as ½AR 0 ¼ 0 in (3b), we have where Substituting in ½R ¼ R tot À ½AR and solving for ½AR gives the usual expression (see also [34], for example) for the concentration of ligand bound at equilibrium as Taking two ligand concentrations, ½A 1 and ½A 2 , and the corresponding output measurements, ½AR 1 and ½AR 2 gives These equations contain the two unknown parameters R tot and K A ¼ k aþ =k aÀ .Solving for these gives a unique solution and hence these parameters are identifiable from a single dose-response curve.In fact, only two points on the curve are needed, theoretically.Once these are known, we conclude from (77) that only a single parameter, either k aþ or k aÀ , remains to be found.This can be obtained from time course data.Using only one of the parameter combinations in (75), we find Hence, we conclude that, using equilibrium data (a doseresponse curve) together with a single set of time course data, it is possible to identify all three model parameters ðk aþ ; k aÀ ; R tot Þ.

Washout experiments
We consider using washout experiment data to identify dissociation parameters.In a washout experiment, the ligand is removed from the system (usually once equilibrium has been reached), hence, we set ½A ¼ 0 in the model given in equations (3a-d).This gives As the concentration of free receptor is unknown at the start point of washout, we write the initial conditions as where the w refers to the value when washout begins, at time t ¼ 0. The output remains unchanged The unknown parameters in this model are k aÀ and R tot .
We note that, it is clearly possible to solve the ODE for ½AR, specifically giving and use the result to determine identifiability directly.However, we refrain from this here and continue with our SIA methodology applied to the ODE system in order to highlight the general process.
Here, the Taylor series approach is straightforward (as outlined in Sect.2.2); no reformulation of the ODE system is required and the unknown initial conditions are naturally incorporated into the analysis.Since the number of variables n ¼ 2, we need to determine a maximum of three Taylor coefficients (see (15)).Calculating the first of these coefficients gives which clearly gives no information about k aÀ , but does provide the value of ½AR w .The second coefficient is and so we find k aÀ to be identifiable.The final coefficient is calculated as which gives no further information, hence, only k aÀ is identifiable from washout data (as we would expect, given (83)).Since R tot does not appear in (84)-( 86), it is clearly not identifiable from the washout model alone.
Combining the newly established identifiability of k aÀ with the results we obtained from a binding timecourse, namely the parameter combinations in (75), we find that the remaining parameters, k aþ and R tot , are now identifiable.
Hence the system, considering the combination of both experiments, is now globally identifiable, using experimental data for association and washout for a single ligand concentration.

Multiple time courses
Next, we consider the case where, instead of one time course, we have multiple sets of time course data, each with a different ligand concentration.Each of these will individually give the identifiable parameters, as stated in equation (75) for their corresponding concentration of ½A.
That is, the identifiable parameter combinations are given by for i ¼ 1; 2; :::, for the number of time courses being considered.As fitting may be performed on all sets simultaneously, we analyse the corresponding system as a single system.For example, for two time courses we have Setting fðpÞ ¼ fðe pÞ results in all parameters being successfully identified.This is shown easily by considering the following system: From (89b,d), we see that and so k aþ is identifiable.Then (89a) gives R tot ¼ g R tot and (89b) gives k aÀ ¼ g k aÀ , so that all three parameters are identifiable.We conclude that the model is fully identifiable from just two time courses.

Competition binding model
Next, we consider the model for a monomeric receptor binding with two ligands in a competition binding scenario, with a labelled ligand A and an unlabelled ligand B. This is described by the schematic and the related system of equations is given in (36a-e).In Sect.3.1, we performed the identifiability analysis considering the parameters p ¼ ðk aþ ; k aÀ ; k bþ ; k bÀ ; R tot Þ, and concluded that from a single time course the parameter k bÀ is uniquely identifiable, as well as the parameter combinations A simple analysis towards establishing identifiability of all parameters is suggested by the scenario discussed in [44] where k aþ and k aÀ are already known from other experiments.For example, we can consider ½B ¼ 0, whereby there is no competition, and use the monomeric receptor model of Sects. 2 and 4.1 to ensure identifiability of these parameters.Then we may treat k aþ as known in the first row of (90), meaning that R tot is identifiable.Treating k aþ and k aÀ as known in the second row of (90), we see that k bþ also becomes identifiable.Hence, with prior knowledge of k aþ and k aÀ , the original system is globally identifiable from a single timecourse with ½B 6 ¼ 0. Continuing with our detailed tutorial approach, we now consider experiments which may be used to establish identifiability of k aÀ , k aþ , k bþ and R tot without the prior knowledge of the binding parameters for ligand A. Again we consider equilibrium concentration-response, a washout timecourse curve or a second timecourse.

Saturation curves
Here, we combine the above results with an equilibrium concentration-response curve.It follows from system (36a-e) that at equilibrium, we have the relations where K A ¼ k aþ =k aÀ and K B ¼ k bþ =k bÀ .Moreover, we have the conservation law Combining these and solving for ½AR gives the equilibrium concentration of the measured ligand bound as Assuming we have a concentration-response curve, with a fixed concentration of ligand ½B, then we take two points on this curve, giving This results in a system of two equations for K A and R tot , which has the following solution: Considering further points on the dose-response curve gives no extra information.Hence from a single concentration-response curve the only identifiable parameter is R tot , since the expression for K A still depends on the unknown parameter K B .We combine the above with the results from a time course dataset, as given by the parameter combinations in (90).Since R tot and k bÀ are now known, the unknowns are given by p ¼ ðk aþ ; k aÀ ; k bþ Þ.We consider fðpÞ ¼ fðe pÞ, where f is given in (90).The first equation can be solved as thus k aþ is now identifiable.The other two resultant equations, using that k bÀ and k aþ are known, are There are two possible solutions of these two equations for k aÀ ; k bþ , given by Hence the system is only locally structurally identifiable, but not globally.This could be circumvented with prior knowledge about the parameters.If we have prior knowledge that k bÀ [ k aÀ , then the first solution in (99) cannot be satisfied since all parameters must be positive.This implies that the only solution that fits within the requirements is the second solution, and hence, that all parameters are identifiable.Remarkably, it is shown in [44]  We may also consider an equilibrium concentration-response curve for varying concentration of the competition ligand B, as in [44], in addition to the single time course and the concentration-response for [A] described above.For example, taking three equilibrium experiments with concentrations ([A], [B]) = ð½A 1 ; ½B 1 Þ, ð½A 2 ; ½B 1 Þ and ð½A 1 ; ½B 2 Þ and corresponding readouts / 1 , / 2 and / 3 gives the following system for the equilibrium parameters K A ; K B and R tot : Solving these for the equilibrium parameters, K A , K B and R tot , gives the unique solution and thus determines all three of these parameters to be identifiable.Combining these with time course results, as in (90), we find the model to be globally structurally identifiable.

Washout experiments
In this section, we combine the results from an association time course experiment, in (90), with a washout experiment.To determine identifiability using washout experiment data we set ½A ¼ 0 in equation (36a-c) to simulate the washout of ligand A (a similar analysis is possible for washout of ligand B, or of both ligands).This results in the reduced model For this system, the initial conditions are unknown so we assume The only measured quantity is ½AR and so the output is Assuming that k bÀ is known, the washout model has unknown parameters p ¼ ðk aÀ ; k bþ ; R tot Þ.We again use the Taylor series method.The first Taylor coefficient is given by and the quantity ½AR w is known.The second coefficient is determined by giving that k aÀ is identifiable.Further coefficients of the Taylor series give no new information, and so only k aÀ is identifiable from dissociation data (again, as in Sect.4.1.2,R tot is not identifiable from washout alone).So use of this washout experiment in conjunction with the binding experiment of Sect.3.1 allows us to consider k aÀ as known in (90).We then determine identifiability using the combination of experiments by considering the possible solution, for ðk aþ ; k bþ ; R tot Þ, of the system It is straightforward to show that the system has a unique solution We conclude that the Motulsky-Mahan problem, considering the combination of both one binding timecourse and one washout timecourse, is now globally identifiable.

Multiple time courses
There are two ways in which we can use multiple sets of time course data to determine identifiability of this model.Taking the coefficients as stated in equation (90) we consider multiple time courses with either several concentrations of A or several concentrations of B, giving either ð107Þ for i ¼ 1; 2; ::: .If we study the first case, having time courses for two A concentrations, we find the identifiable parameter combinations, in addition to the single identifiable parameter k bÀ , are given by : Solving fðpÞ ¼ fðe pÞ, we find that k aþ and R tot are identifiable and that k aÀ and k bþ again satisfy (99).The identifiability properties of the Motulsky-Mahan system combining timecourses for two different values of [A] are the same as those for the single timecourse plus the concentration-response curve for [A].When we instead consider having time courses for two concentrations of B, we obtain the vector of identifiable parameter combinations (in addition to k bÀ , which we know is identifiable from a single timecourse) Solving fðpÞ ¼ fðe pÞ this time gives that all parameters are identifiable from these two time courses.

Pre-dimerised G protein-coupled receptor binding
In Sect.3.2 we explored identifiability for a model of dimeric receptor binding with a single ligand.The schematic for this is given by and the system of equations is given in (61a-e).We found that from a single time course we have no identifiable parameters.Recall that this analysis gave the following identifiable parameter combinations (see (66))

Saturation curves
We first consider the combination of (110) with information from an equilibrium concentration-response curve which is parameterised by equilibrium parameters, The expression for the concentration of ligand bound at equilibrium, A b was determined in our previous work [55], and is given by Taking three points on the dose-response curve, that is, three different concentrations of ½A i , we have where ½A b i denotes the corresponding measurement for the concentration of ligand ½A i .Equation (112), with i ¼ 1; 2; 3, is a system of three equations which can be solved for K A , a and R tot .This can be done by using a symbolic equation solver (for example, in MATLAB [40] or Mathematica [58]), and using MATLAB Symbolic Toolbox, we find a unique solution to (112), using just three points.The expressions are extremely lengthy and impractical to write down, and therefore we refrain from doing so here.We now treat K A , a and R tot as known quantities, together with identifiable parameter combinations (110) to determine identifiability of the parameters in p ¼ ða þ ; a À ; k aþ ; k aÀ Þ.The first of the equations, after setting fðpÞ ¼ fðe pÞ in (110), yields Clearly, k aþ is identifiable.Combining this with the known equilibrium parameter K A ¼ k aþ =k aÀ , we find that k aÀ is also identifiable.The second and third equations from fðpÞ ¼ fðe pÞ can be simplified to , and thus all four parameters a þ ; a À ; k aþ and k aÀ (in addition to R tot , from the equilibrium analysis) are identifiable.

Washout experiments
Here we consider washout experimental data.The corresponding model for washout of the ligand is given by setting ½A ¼ 0 in equations (61a-e), giving The initial conditions are unknown so we assume where subscript w denotes the value at the start of washout, and ½AR w and ½ARA w are unknown.The output remains as We again use the Taylor series method to determine identifiability in this section.Through repeated differentiation of y and substitution of the initial conditions we obtain the vector of coefficients as Since R tot does not appear in f 1 here, we note that it is not identifiable from the washout model alone, as in Sect.4.1.2.Further, while the parameters a À and k aÀ are sought, we have introduced new parameters ½AR w and ½ARA w , the initial conditions, which are also unknown and are intertwined in the identifiable parameter combinations in (117).Although we do not require them to be identifiable, the analysis requires us to consider them.In Sects.4.1.2and 4.2.2, the analysis was simpler given that the new parameter ½AR w was immediately identifiable as the measured output.At this point, we may define p ¼ ða À ; k aÀ ; ½AR w ; ½ARA w Þ for the washout experiment and proceed by attempting to solve f 1 ðpÞ ¼ f 1 ðe pÞ to determine the identifiability of individual parameters.Thereafter, we could return to the association timecourse result (110) to determine identifiability of those individual parameters.Given the level of complexity seen in ( 117), we would use symbolic computation here.Alternatively, our computation could consider p ¼ ða þ ; a À ; k aþ ; k aÀ ; R tot ; ½AR w ; ½ARA w Þ for the two experiments (association and washout) combined.Then the combined identifiable are groupings given by Now when solving f 2 ðpÞ ¼ f 2 ðe pÞ, we find (by symbolic computation, see Appendix 4) the unique solution Hence we conclude that the combination of association (61a-e) and washout (116a-e) results in all five of the parameters ða þ ; a À ; k aþ ; k aÀ ; R tot Þ being globally identifiable.

Multiple experiments
Next, we consider association timecourse data obtained from experiments each with a different ligand concentration [A].We aim to determine the minimum number of concentrations needed to ensure full identifiability.Each of these experiments yields identifiable parameter combinations, as obtained in (110), with ½A i as the concentration for experiment y i .This gives the identifiable parameter combinations Upon assuming data for two experiments, we choose i ¼ 1; 2 to give the vector of identifiable combinations as To determine identifiability we again set fðpÞ ¼ fðe pÞ and solve for p ¼ ða þ ; a À ; k aþ ; k aÀ ; R tot Þ.It is a matter of simple algebraic manipulation to find that So we conclude that only two data sets are required to ensure that the system is globally structurally identifiable.

Discussion
Structural identifiability analysis (SIA) is an often-overlooked element of modelling of biological systems [12].The notion of identifiability is well known and appreciated, but in practice the complexity of the calculations that are required to draw conclusions regarding the identifiability of a given ODE system is often a barrier.While the ''classical'' SIA methods of transfer function, Taylor Series and similarity transformation have been applied to a number of pharmacokinetics models in the literature, SIA is largely absent from receptor theory and analytical pharmacology studies.Here, we have introduced SIA methodology to receptor theory via application of these three classical methods to three widely adopted ligand-receptor binding schematics of biological importance.Our analysis has yielded new identifiability results for single-timecourse receptor theory outputs, plus a significant and crucial focus on approaches to mitigating non-identifiability via the addition of further experiments.In addition, the article provides a pedagogical, tutorial-style introduction to formal identifiability analysis, aligned with the aim of bringing SIA to a broader audience [7].
Our key results include a formal SIA verification that, for the model of ligand A binding monomeric receptor, the quantities k aþ R tot and k obs ¼ k aþ ½A þ k aÀ (the so-called observed on-rate [34]) are globally identifiable, see expression (6).For this monomeric receptor model, we also confirm mathematically the intuitively known fact that all kinetic parameters and the receptor concentration become identifiable when adding a washout experiment, an equilibrium saturation binding model or simply when using timecourses for two ligand concentrations.
For the Motulsky-Mahan competition binding model (36a-e), it was already known that if the total receptor R tot and labelled ligand constants k aþ ; k aÀ are known, then the unlabelled ligand constants k bþ and k bÀ are theoretically identifiable from a single [AR] timecourse [15,44].Our SIA (see expressions (43,44)) shows that k bÀ is in fact identifiable without the need for a priori knowledge of any constants.Furthermore, if k aþ is known, then R tot is also identifiable without knowledge of any other parameters, or if R tot is known, then k aþ is identifiable.From these results, it is clear that a single timecourse may yield more practical, quantitative parametric information than previously thought.In addition, we have shown that SIA enables a formal strategy for constructing an identifiable system when also considering washout experiments and/or multiple ligand concentrations (Sect.4.2).Recent computational studies of the Motulsky-Mahan model have focused on questions of practical identifiability and parameter estimation [18,49] and have noted a relationship between binding timescales and estimation reliability.Further investigation into this relationship will benefit from the analytical results and methods presented in the current work.
We have also shown new results from SIA applied to a model of GPCR dimers in Sect.3.2.This model has been previously used with experimental data for total bound ligand to partially quantify the important effect of cooperativity (reporting an equilibrium parameter) across a dimer [41] without discussion of parameter identifiability properties.Our new analysis indicates that no model parameters are identifiable from a single binding timecourse.However, when using multiple ligand concentrations, all kinetic parameters and the total receptor concentration are identifiable.These results that were so far unknown provide both practical guides for the estimation of kinetic cooperativity and an extension of the recent theoretical study of cooperativity and dimer binding dynamics given in [55].
The models we have considered have been low-dimensional (at most third-order) and linear, which is typical of many ligand-binding models in receptor theory.For such models, the implementation of the three classical SIA methods is tractable.Given the relative conceptual simplicity of these approaches compared to more recent methods developed for larger biology and systems biology models [3,12,46,48], the introduction of SIA to receptor theory via these three methods has been shown to be viable, although we remark that the Taylor Series approach is cumbersome in some cases, benefiting from symbolic computation tools.The transfer function method is relatively straightforward, and reflects earlier use of the Laplace Transform in textbook PK parameter estimation discussion [19].The Taylor Series approach and a modified similarity transformation method suitable for nonlinear models [10,50] are potentially suitable for second-order nonlinear binding models such as those arising in ligandinduced receptor dimerisation [56], and simple nonlinear models for receptor-mediated cell responses via kinetic operational models of agonism [26].To bridge the gap to higher-dimensional models of interest in receptor theory (including binding of allosteric modulators [24], more detailed operational models [26] and G protein activation [59]), recent novel algorithmic and computational approaches including Exact Arithmetic Rank [46], inputoutput method [3,31] and singular value decomposition of sensitivity matrices [48] appear to be promising methods.
We conclude by proposing the following studies as future work.
1. Apply SIA to a binding model for ligand-induced receptor dimerisation.Recent analysis of this nonlinear ODE model has shed new light on dynamic cooperativity effects across the dimers, and the model has been validated by fitting it to real timecourse data [56].SIA is required to determine the theoretical identifiability of the model parameters (in preparation [57]).2. Apply an identifiability analysis to the Motulsky-Mahan model combining both the structural identifiability results presented in the current work and the recent practical identifiability and estimation results in [15,18,49] to derive an overall strategy for informing parameter estimation studies for this widely-used model.3. Perform a bridging-the-gap analysis which implements the Exact Arithmetic Rank [46], input-output method [3,31] and singular value decomposition method [48] to the models in the current work, the ligand-induced dimerisation model and kinetic operational models, to compare ease of implementation and computational cost.
the system to be controllable and observable.In the full system (36a-e) we have G ¼ ð0; 0; 0Þ T , and so the controllability matrix is In this case, the dimension is n ¼ 2, and the controllability matrix is given by which has rankðCÞ ¼ 2, and so the system is controllable.The observability matrix is found to be which also has rankðOÞ ¼ 2, therefore, both controllability and observability conditions are met.Next, we assume that there exists a transformation matrix and assume that conditions (32a-e) have to hold.Then, step by step, we draw conclusions from these conditions, and find the entries in matrix T. We first note that (32b) is satisfied for any T as x 0 ¼ x0 ¼ ð0; 0Þ T .We now check condition (32e), which becomes Then, condition (32d) implies that 1 0 which leads to The first of these equations, since ½A is known, yields which is the first identifiable parameter combination.
Next, (131) can be written as So far, we have found that Next, for condition (32c) to hold needs to be satisfied.We multiply both sides by À1, and write (135) as where and We now equate the entries of matrices M and N. First we solve M 12 ¼ N 12 for t 21 , giving We equate M 11 ¼ N 11 , which together with the above expression for t 21 , becomes Since we have already determined that g k aþ g R tot ¼ k aþ R tot (see (132)), this results in Hence we find which is the second identifiable combination.Equating M 22 ¼ N 22 gives With g k aþ g R tot ¼ k aþ R tot , this can simplified to confirming that the parameter k bÀ is identifiable.Finally, equating M 21 ¼ N 21 , gives, after some simplification (again using MATLAB Symbolic Toolbox) which gives a fourth identifiable parameter combination.We note that the transformation matrix is given by As detðTÞ 6 ¼ 0, we find that condition (32a) holds.In summary we have k bÀ as identifiable, and identifiable parameter combinations That is, the same parameter combinations are found to be identifiable as for the previous methods (see (44) and ( 59)).
Appendix 2 Taylor series approach to GPCR homodimer model Here, we apply the Taylor Series method to the GPCR dimer model.In (73a-d), we established that for the GPCR homodimer model, a vector of identifiable combinations is then given by where By ad hoc algebraic manipulations, largely using MATLAB Symbolic Toolbox, we find that showing that the combinations in (74) are recoverable from (147a-g) and hence identifiable.A sample MATLAB code is given in Supplementary Materials.
Before we determine identifiability for the system, we first check whether the controllability and observability conditions are satisfied.Using (29) the controllability matrix for this system is found to be and the observability matrix is found, using (30), to be As rankðCÞ ¼ rankðOÞ ¼ 2 we conclude that the system is both controllable and observable, hence we continue with the identifiability analysis.
Assuming the existence of a linear transformation matrix and that conditions (32a-e) hold, we begin to draw conclusions.We first note that xð0Þ ¼ 0 implies that condition (32b) automatically holds.We now apply condition (32e), giving From this we determine We next apply condition (32d), giving 1 À 2t 21 2 À 2t 22 From the bottom row, we find Then the top row gives the first identifiable parameter combination as With (156) the transformation matrix now becomes Applying condition (32c) and determining the left-hand and right-hand side we obtain and We solve M 21 ¼ N 21 for t 22 , which gives We note that, while it is possible to begin with a different matrix entry, and therefore have a different expression for t 22 , this will still give the same identifiability results.We substitute the expression for t 22 into the remaining matrix entries.Solving M 11 ¼ N 11 , we find, after some simplification Equating M 22 ¼ N 22 and M 12 ¼ N 12 , we obtain respectively.Hence, we find the identifiable parameter combinations for this model given in the vector The transformation matrix is given by and so we confirm that condition (32a) holds.Note that the combinations listed in (164) and ( 74) are not identical.However, since the second and fourth entries in (164) sum to give we conclude that the identifiable parameter combinations are the same as those in (66) and (74), namely In Supplementary Materials, we include a MATLAB code which, when run after the code for the Appendix 2 computations, returns a unique solution to the equation f 2 ðpÞ ¼ f 2 ðe pÞ for p ¼ ða þ ; a À ; k aþ ; k aÀ ; R tot ; ½AR w ; ½ARA w Þ and f 2 given by (118).

Fig. 2
Fig. 2 Time courses for observed output a c ðtÞ and non-observed a b ðtÞ for the two-compartment PK model (1a-c), for the two parameter sets given in Table1

Fig. 5
Fig.5Three sets of parameters are used to plot the solution to equations (61a-e).All three parameter sets give the same measured output curve, A bound .However, non-identifiability can be seen in the

k
aþ g k aþ ½A g R tot

Appendix 4
Combined association and washout computation for GPCR homodimer model

Table 3
The values of three different parameter sets that are used to plot Fig.4 Table 4 together with ½A ¼ 10 À8 M that k bÀ [ k aÀ if and only if [AR](t) is monotonic.Therefore, an experimental timecourse readout showing no peak in [AR](t), together with an equilibrium concentration-response curve for ligand A results in identifiability of all five parameters ðk aÀ ; k aþ ; k bÀ ; k bþ ; R tot Þ.For the same combination of experimental data, but with non-monotonic [AR](t), only ðk aþ ; k bÀ ; R tot Þ are identifiable.