Abstract
Background
A doseresponse curve depicts the fraction of bound proteins as a function of unbound ligands. Doseresponse curves are used to measure the cooperativity degree of a ligand binding process. Frequently, the Hill function is used to fit the experimental data. The Hill function is parameterized by the value of the dissociation constant and the Hill coefficient, which describes the cooperativity degree. The use of Hill's model and the Hill function has been heavily criticised in this context, predominantly the assumption that all ligands bind at once, which resulted in further refinements of the model. In this work, the validity of the Hill function has been studied from an entirely different point of view. In the limit of low copy numbers the dynamics of the system becomes noisy. The goal was to asses the validity of the Hill function in this limit, and to see in what ways the effects of the fluctuations change the form of the doseresponse curves.
Results
Doseresponse curves were computed taking into account effects of fluctuations. The effects of fluctuations were described at the lowest order (the second moment of the particle number distribution) by using the previously developed Pair Approach Reaction Noise EStimator (PARNES) method. The stationary state of the system is described by nine equations with nine unknowns. To obtain fluctuationcorrected doseresponse curves the equations have been investigated numerically.
Conclusions
The Hill function cannot describe doseresponse curves in a low particle limit. First, doseresponse curves are not solely parameterized by the dissociation constant and the Hill coefficient. In general, the shape of a doseresponse curve depends on the variables that describe how an experiment (ensemble) is designed. Second, doseresponse curves are multivalued in a rather nontrivial way.
Background
The Hill function is frequently used to infer the degree of cooperativity of the chemical reaction in which ligand molecules bind to a protein [1]. Often, the binding of a ligand increases the association rate for the binding of the next ligand. Such reactions are said to be (positively) cooperative. There are examples of cooperative reactions in cell biology. The classical example is the binding of oxygen molecules by hemoglobin [1]. Other perhaps less wellknown examples would be parts of the Notch signaling and 30 S ribosome assembly processes [2], as well as the assembly of cholesterolsphingomyelin complexes [3]. Also, the noise characteristics of various ligand binding reactions were studied theoretically in [4] and some of the experimental systems could be classified as cooperative reactions. A cooperative reaction builds a final complex successively. If strong cooperativity is present, the dynamics of the system can be studied using Hill's model, at least to a first approximation [5].
Hill's model is a grossly simplified version of reality. The model is constructed by assuming that binding and unbinding of ligands occur in one step as
where C_{0 }denotes a protein that binds ligands A, and C_{h }is the ligandprotein complex. The Hill coefficient h describes the number of binding sites on the protein. Both the forward and the back reactions are allowed.
Strictly speaking, the Hill coefficient in Hill's model (1) is a stoichiometry coefficient and should be an integer number larger than zero. However, in the calculations that follow, h will be allowed noninteger values. Thus in the context of this work the Hill model should be understood more from a model average perspective, where the Hill coefficient is an effective parameter.
An important quantity related to Hill's model is the fraction of the proteins that are bound
In particular, the dependence of φ on the amount of unbound ligand in the system a is of considerable interest, and is referred to as a doseresponse curve. A function frequently used to fit a doseresponse curve is the expression derived by Hill, the socalled Hill function, given by
where c_{0}, c_{h}, a are used to denote the amounts of unbound proteins, bound proteins, and free ligands, respectively. Please note that the Hill function is only parameterized by K and h. When fitting experimental data to extract K and h, it is useful to allow h to be a real number. Also, the Hill function is used frequently in theoretical studies to model cooperativity effects.
In general, c_{0}, c_{h }and a can denote average particle numbers, particle concentrations or partial pressures. It really depends on the types of experiments one wishes to describe. The dissociation constant is essentially controlled by the ratio of the forward and the backward reaction rates.
The original Hill's model is unrealistic since a truly multiparticle reaction with a high Hill's coefficient would be a very unlikely reaction event. The probability that all required ligand molecules meet at the right place, at the right time, is very small. The model was already criticised by Hill himself [6,7]. Subsequently, more realistic models were suggested in a series of studies: Adair [8]; Monod, Wyman, Changeux [9]; and Koshland, Nemethy, Filmer [10]. The difference between the models was critically investigated on the mean field level in [5], which confirmed Hill's original claim that the Hill equation can be used in a case of strong cooperativity when intermediate states are shortlived. For a reaction set that appears strongly cooperative as in (1), the Hill coefficient provides a rough measure of the cooperativity degree of the reaction.
Despite the problems discussed above, the use of Hill's model has some merits [1], and the Hill equation is used frequently in many fields as discussed in review article [11]. Accordingly, in this work, Hill's model will be taken as a basic standard for describing multiparticle (cooperative) reactions. The validity of the model has been extensively investigated previously. The conditions for safe usage of Hill's model can be easily verified.
From now on, it will be assumed that the Hill model under investigation is a valid alias for a more complicated multiparticlelike reaction scheme. The focus will be on investigating the correctness of the resulting Hill's function φ_{H}(a) in a low particle number limit. The ultimate goal of this study is to investigate in what ways the effects of the noise related to the low copy numbers affect the form of the doseresponse curve predicted by Hill. Please note that such a goal enforces consideration of a closed system. For an open system, where injection and the decay of particles are allowed, one cannot use the Hill function at all.
Results and discussion
Model description
The fundamental quantity we wish to understand is the fraction of bound proteins φ in a situation when particle numbers are low. This is done by considering a closed system in a well mixed regime. In such a situation it is sufficient to count the particles. In the following, n_{0}, n_{h}, and n_{A }will denote the number of C_{0}, C_{h}, and A particles respectively. A stochastic model will be considered with the forward reaction rate α and the back reaction rate β. The rates have the dimension of inverse time. Owing to the stochastic nature of the model, the particle numbers will fluctuate. The ensemble averages of fluctuating quantities will be denoted by ⟨.⟩. Accordingly, particle amounts will be expressed in terms of average particle numbers, c_{0 }= ⟨n_{0}⟩, c_{h }= ⟨n_{h}⟩, and a = ⟨n_{A}⟩. In such a case the dissociation constant in equation (3) is precisely given by
The expression for K in (4) can be obtained from the stationary state equations that describe the system in the mean field limit. Use of equations (2729) and (30) in the methods section leads to the desired result. Strictly speaking, the variable K is not a dissociation constant, but it can be related to it by trivial rescaling by the volume of the system.
For any type of initial conditions the dynamical system at hand will reach equilibrium. The focus will be on investigating the equilibrium state of the model, which in turn will enable us to compute the dose response curve φ (a).
Analytical description of system is possible
The central technical result of this paper is the derivation of the nine (nonlinear) equations (513) with nine unknowns. These equations describe the equilibrium state of the model. The derivation of the equations is described in the methods section. The equations can help in analytical understanding of the problem.
The first three stationary state equations are given by
In equation (5), and in the following, the symbol χ with a subscript denotes a correlation function. Correlation functions were introduced previously (Konkoli, Z.: Multiparticle reaction noise characteristics, submitted) and describe fluctuations. The situation when all χ = 1 corresponds to the mean field limit, where the effects of fluctuations are absent. It is easy to see that in such a case equations (57) combine to give the classical Hill function in (3). However, the correlation functions do not equal one in general, and the expression for the Hill function in equation (3) might be invalid.
Equations (6) and (7) express the fact that the total number of protein complexes (with and without ligands) P_{0}, and the total number of ligands in the system (both free and bound) L_{0}, cannot change over time. Averages ⟨P_{0}⟩ and ⟨L_{0}⟩ need to be used; depending on an ensemble, these quantities might be stochastic. It ultimately depends on how the system is prepared during an experiment.
The remaining six equations feature correlation functions heavily. The first three are
and are obtained from analysis of the dynamics that brings the systems to a stationary state. The last three equations are the conservation laws that express the fact that initial fluctuations in P_{0 }and L_{0 }cannot change over time:
The nine equations with the nine unknowns (513) are the central result of the paper.
The equations are nonlinear and fully describe the stationary state of the system
when the effects of particle number fluctuations are taken into account. The observables
of interest (average numbers of particles and correlation functions) are implicit
functions of the ensemble properties ⟨P_{0}⟩, ⟨L_{0}⟩,
The equations are not exact. They were derived using the Pair Approach Reaction Noise Estimator (PARNES) method introduced previously (Konkoli, Z.: Multiparticle reaction noise characteristics, submitted). The PARNES method works by approximating higher order moments of a particle number distribution by second order moments. Should the need arise, the method can be easily extended beyond the pair approach level.
The PARNES method is based on the usage of correlation forms. The correlation forms are used in studies of spatially extended diffusion controlled reactions [12]. They are employed to close the hierarchy of manypoint density functions. In the present work, the particular methods discussed in [13] were adopted to study a well mixed reaction system. Because a second quantization formalism is used, the PARNES approximation is naturally expressed as a closure relationship for factorial moments of a particle number distribution. The implementation of the closure procedure is shown in the methods section. There are other ways to perform the closure [4,1418].
Clearly, once moments are given it should be possible to work backwards and extract the form of the particle number distribution function. This is a rather nontrivial problem and will be studied elsewhere. Essentially, the PARNES approximation is an expansion around the Poisson distribution. For χ ≈ 1 the distribution function is Poissonlike. Situations with χ < 1 and χ > 1 describe sub and supraPoisson regimes respectively.
The Hill equation is valid for large copy numbers
It is possible to see that when particle numbers become large the correlation functions
approach the mean field limit in which all correlation functions are equal to one.
For example, by neglecting the ah^{2}c_{h}, ⟨P_{0}⟩ and hc_{h }terms in equations (11), (12) and (13) respectively, and assuming that
A danger of inferring an incorrect Hill's coefficient
The issue is whether all solutions of the central equation system are such that φ can be expressed solely as a function of a. If this is the case then there is only one equation to use, and there should be no ambiguity regarding the proper choice of Hill's coefficient. By inspecting the form of the central equations it can be seen that this is not the case in general. For example, depending on the procedure used to compute the points in the plot that depicts φ (a), many curves can be obtained. Equivalently, in more technical terms, for a given reaction system, repeating the experiment to determine φ (a) with different ensemble setups (the ways the system is prepared), one can obtain different curves for φ (a). Fitting the curves to φ_{H}(a) would result in different Hill's coefficient for each curve. Thus, the fact that the central equations depend on ensemble properties has far reaching consequences when it comes to extracting the correct Hill coefficient from experiments.
Numerical tests
The question is how much the effects of noise affect the shape of doseresponse curves. To address this question the nine equations were solved numerically for relatively low copy numbers of the protein that binds ligands. Figures 1 and 2 shown that φ is not solely a function of a, but depends on the characteristics of the ensemble as suggested. The figures describe the Poisson and pure ensembles respectively. The curves in the figures clearly depend on the way that is used to prepare the initial state of the system.
Figure 1. Fraction of bound proteins (Poisson initial state). A doseresponse curve (the fraction of the bound proteins φ plotted as a function of a) for a Poissonlike ensemble:
Figure 2. Fraction of bound proteins (pure initial state). Does response curves for the system prepared in a pure state:
Analysis of both figures shows that for large particle numbers the mean field result (the Hill function) is obtained. This is expected, since the mean field description should be correct for large copy numbers. However, in general, the discrepancy from the mean field case can be significant. For Poissonlike initial conditions the reference curve is approached from below. In the case of pure initial states, the reference curve is approached from above (below) for high (low) values of a.
For pure initial states, and in the intermediate regions of a, φ curves are much steeper that the corresponding Hill function. Please note that the curves for pure states are multivalued since for a given value of a there can be more than one value of φ (e.g. all thin curves in Figure 2 for values of a slightly greater than one are multivalued). Similar behaviour is observed for Poissonlike initial states but the onset occurs at smaller values of a (e.g. the dotted line in Figure 1). The question is whether such behavior is an artefact of using the PARNES approximation.
Figure 3 depicts φ (a) obtained by an exact diagonalisation of the master equation. The figure shows that φ (a) is indeed multivalued. The exact solutions exhibit richer behavior than is predicted by the PARNES method. It is very likely that the erratic alternation of points has to do with the fact that not all ligands can be fully absorbed by the receptors. For example, assume that one observes a snapshot of the system dynamics where all proteins in the system have bound all ligands. If one adds more ligands to the system, any number in range from 1 to h  1, exactly that number of ligands will never be bound by the receptor proteins. A similar effect was observed in a related study [19]. Such effects cannot be explained directly by usage of the PARNES method. The PARNES method can describe such behavior only qualitatively.
Figure 3. Fraction of bound proteins (pure initial state), exact result. Exact dose response curves for a system in pure states. As in Fig. 2 the thickest full line is the reference Hill curve. Thinner curves were generated by direct diagonalisation of the master equation. The thinner full lines are obtained for fixed value of P_{0 }and looping values of L_{0}. For each point (L_{0}, P_{0}) the master equation was solved numerically and observables of interest were computed. The full line is for P_{0 }= 2. The dashed line is obtained for a much larger number of receptors P_{0 }= 8. This figure shows that exact dose response curves are multivalued. Since not all points are physical, the points were connected using linear interpolation to guide the eye. The dose response curves obtained in such a way are rather erratic. Furthermore, the multivalue character is not an artefact of using linear interpolation. There are many physical points with nearly identical values for a having many distinct values for φ.
Figure 4 depicts φ as a function of L_{0 }for a pure ensemble. From a theoretical point of view the dependence of φ on a is of interest, but φ is more likely to be plotted as a function of L_{0 }in experimental work. Please note that φ (L_{0}) is a single valued function. However, the curve depicting the exact dependence of φ on L_{0 }is not smooth. The notion of the curve is to be understood by interpolating between allowed points since only integer values for L_{0 }make sense for a pure ensemble. The curve obtained by using the PARNES approximation follows the exact result much more closely than the mean field curve.
Figure 4. Fraction of bound proteins; L_{0 }dependence. The fraction of the bound proteins φ is plotted as the function of free ligands in the system L_{0 }for the pure state. All curves were obtained for P_{0 }= 2. The thickest full line is the mean field result. The thinner full line is obtained using the PARNES method. The dashed curve is obtained by exact diagonalisation of the master equation. Please note that the PARNES curve (thin full line) agrees best with the exact result (dashed line).
Conclusions
Many dangers have already been recognized in using the Hill function to fit experimental data. The difficulties discussed so far in the literature are mostly related to the fact that the Hill model is only an approximation of a more complicated reaction scheme. This work points to a yet another danger, but in terms of principles.
The findings of this work point to the fact that one should be careful in using the Hill function to fit experimental data when the number of particles in the system is low. The actual dependence of φ on a is much more complex than predicted by the Hill function φ_{H}(a). First, doseresponse curves depend on the way the experiment is done. Repeating the experiment with different ensemble properties could result in a number of distinct curves. Accordingly, equally many values for the Hill coefficient could be extracted. Second, doseresponse curves are multivalued in a rather nontrivial way, which has to do with the fact that some ligands will always be unbound, depending on the number of ligands in the system.
The discrepancy between fluctuationcorrected doseresponse curves and the Hill function has nothing to with a fundamental flaw in the Hill model itself. The features are rather generic. Similar behaviour is likely to be observed for any more realistic model of ligand binding.
The nine equations obtained in this work could aid experimental studies in which the Hill coefficient is measured. Clearly, to obtain the correct value for the Hill coefficient, one needs to use the correct curve. The nine equations that define doseresponse curves could be investigated further to obtain analytical approximations for fluctuationcorrected doseresponse curves.
This work can be extended in many ways. The uniqueness conditions for the equations have not been investigated yet. Preliminary numerical investigations show that the structure of the solutions is rather complex, since Mathematica solver had to be finetuned to find the solutions. Also, the nine equations allow for nonphysical solutions with negative densities or negative correlation functions. This problem can be solved by proper parameterization of the densities. The question is whether some of the features observed here are an artefact of the "all or none" reaction principle that is intrinsic to Hill's model. For example, it is not clear whether the multivalue character of dose response curves will still be observed in more realistic ligand binding models. Some of the issues discussed above will be investigated in forthcoming publications.
Methods
Mapping to quantum field theory
The problem at hand is stochastic and can be described by a master equation:
where ∂_{t }denotes the time derivative, and c = (n_{0}, n_{h}, n_{A}) is a configuration of the system specified by the number of free proteins, ligand proteincomplexes and free ligands. The states c[+,, +] and c[,+,] are defined by
where any combination of the plus and the minus signs can be picked at will. The particle number probability distribution function P(c, t) defines the probability that the system is found in a configuration c at a time t. Please note that the equation contains binomial coefficients that count ways of choosing clusters of h particles.
The quantities of interest are observables of the type
where f is an arbitrary function of state c. In principle, to compute the averages using (16) is hard. Such a procedure would require the direct solution of the master equation, which is computationally rather demanding. To avoid using equation (16), the equations of motion for the observables of interest will derived. Once in place, these equations of motion can be studied directly. To derive the equations, the problem is mapped on to a quantum field theory using the standard techniques [20]. Thereafter, it is possible to derive the desired equations of motion in a straightforward manner. Please note that any other approach can be used to derive the equations. The filed theory is used in here since it is a useful bookkeeping device.
The field theory for the problem is constructed as follows. The particle number probability distribution function is used to construct the generating function
where
and the operators in parentheses denote the creation operators for C_{0}, C_{h }and A particles:
The field theory that describes the problem is defined through the expression for the Hamiltonian operator that describes the dynamics:
The requirement for equivalence between equations (14) and (19) fixes the form of the Hamiltonian operator, which turns out to be
Using quantum field theory formalism, the observable in (16) can be calculated as
where the right hand side of equation (21) is evaluated using the standard commutator rules for the operators
and the fact that
Equations of motion
An equation of motion for the observable
The equation follows from (19) and the fact that ⟨1
Using equation (26) with
where
Please note that the equations contain the expression ⟨ĉ_{0}â^{h}⟩, so it appears that we need an equation for that quantity as well. This will be dealt with later.
The fluctuations in the numbers of particles will be described by the second moments of the particle number distribution for all pairs. The equations for the second moments are given by
Conservation laws
The system is closed and five conservation laws can be extracted from the equations of motion. This can be done by taking the appropriate linear combinations of the equations so that the time derivatives vanish. The first two conservation laws are given by
and express the fact that the total number of protein complexes (with and without ligands) P_{0}, and the total number of ligands in the system (both free and bound) L_{0}, cannot change over time. For example, the first conservation law can be obtained by adding equations (27) and (28).
Related to the two conservation laws discussed above it is possible to derive the three additional laws that describe the conservation of fluctuations in P_{0 }and L_{0}:
Please note that the conservation laws involve only quantities that describe the ensemble
that was used to prepare the system. The ensemble is defined by five independent parameters ⟨P_{0}⟩, ⟨L_{0}⟩,
Stationary state equations
The Hill function describes stationary states. Accordingly, the equations of motion will be studied in the long time limit. Requiring that all time derivatives in equations (2729) and (3136) vanish gives the set of four equations
The equations involve expressions for which additional equations of motion need to be derived. Unfortunately, such a procedure results in an infinite hierarchy of equations. To cut the hierarchy, the PARNES approximation is discussed. In technical terms, all expressions that involve a product of three of more operators are approximated by products of the pair correlation functions. The pair correlation functions are defined as
In the strict mathematical sense the PARNES approximation can be expressed as follows
where x, y and z are integers greater than or equal to zero. The accuracy of the PARNES approximation has been investigated on a similar model where it was confirmed that it provides a semiquantitative description (Konkoli, Z.: Multiparticle reaction noise characteristics, submitted). For large particle numbers it is rather accurate. A similar investigation for Hill's model (Figure 4) leads to the same conclusions. The PARNES approximation provides a qualitative description of the stationary state of Hill's model in the low particle number limit.
Finally, using the PARNES method (52), an approximative form of equations (4245) can be derived. Carrying out the procedure, and combining the result with the conservation laws (3741), results in the nine equations with nine unknowns listed in (513), which were introduced in the results section.
The central equations (513) can be obtained roughly as follows. Equation (5) results from the stationary state condition (42), and equations (67) are the first two conservation laws (37) and (38) expressed in a new notation. Equations (810) result from the stationary state conditions (4345). Equations (1113) are derived from the conservation laws for second moments (3941).
Numerical recipe
In the general case, the equations are rather involved and cannot be solved analytically. The numerical procedure for solving the equations naturally suggests itself as follows. First, one solves equations (57) assuming that all correlation functions are one. This gives the first guess for the average particle numbers c_{0}, c_{h }and a. The values obtained are inserted into (813) to evaluate a guess for the correlation functions. The resulting values can be again used again in (57) to obtain even better values for the average particle numbers. The procedure continues until results converge to the fixed point values.
However, the procedure discussed above is numerically unstable in the low particle number limit. The plots in the work were generated by a method similar to the analytic continuation. The equations were solved in the large particle number limit by the method outlined in the previous paragraph, after which the desired point in the ensemble parameter space can be approached incrementally along a line. In every step, the solution from the previous point is used as a guess for the point that follows.
Competing interests
The author declares that they have no competing interests.
Acknowledgements
The financial support of the Chalmers Biocenter and an internal MC2 grant for strategic development is greatly acknowledged.
References

Ferrell JE Jr: Questions and Answers: Cooperativity.
J Biol 2009, 8:157. BioMed Central Full Text

Williamson JR: Cooperativity in macromolecular assembly.
Nature Chem Biol 2008, 4:458465. Publisher Full Text

Radhakrishnan A, Li XM, Brown RE, McConnell HM: Stoichiometry of cholesterolsphingomyelin condensed complexes in monolayers.
Biochim Biophys Acta Biomembranes 2001, 1511:1. Publisher Full Text

Gurevich KG, Agutter PS, Wheatley DN: Stochastic description of the ligandreceptor interaction of biologically active substances at extremely low doses.
Cell Signal 2003, 15:447453. PubMed Abstract  Publisher Full Text

Weiss JN: The Hill equation revisited: uses and misuses.
FASEB J 1997, 11:835841. PubMed Abstract  Publisher Full Text

Hill AV: The cobinations of haemoglobin with oxygen and with caron monoxide.

Hill AV: The Combinations of Haemoglobin with Oxygen and with Carbon Monoxide. I.
Biochem J 1913, 7:471480. PubMed Abstract  PubMed Central Full Text

Adair G: The Hemoglobin System. VI. The Oxygen Dissociation Curve of Hemoglobin.

Monod J, Wyman J, Changeux JP: On nature of allosteric transiitons  a plausible model.
J Mol Biol 1965, 12:88118. PubMed Abstract  Publisher Full Text

Koshland DE, Nemethy G, Filmer D: Comparison of experimental binding data and theoretical models in proteins containing subunits.
Biochemistry 1966, 5:365382. PubMed Abstract  Publisher Full Text

Goutelle S, Maurin M, Rougier F, Barbaut X, Bourguignon L, Ducher M, Maire P: The Hill equation: a review of its capabilities in pharmacological modelling.
Fundam Clin Pharmacol 2008, 22:633648. PubMed Abstract  Publisher Full Text

Kotomin E, Kuzovkov V: Modern aspects of diffusioncontrolled reactions: cooperative phenomena in bimolecular processes, Volume 34 of Comprehensive chemical kinetics. Amsterdam: Elsevier; 1996.

Konkoli Z: Application of Bogolyubov's theory of weakly nonideal Bose gases to the A+A, A+B, B+B reactiondiffusion system.
Phys Rev E 2004, 69:011106. Publisher Full Text

Elf J, Ehrenberg M: Fast evaluation of fluctuations in biochemical networks with the linear noise approximation.
Genome Res 2003, 13:24752484. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

GomezUribe CA, Verghese GC: Mass fluctuation kinetics: Capturing stochastic effects in systems of chemical reactions through coupled meanvariance computations.
J Chem Phys 2007, 126:024109. PubMed Abstract  Publisher Full Text

Lee CH, Kim KH, Kim P: A moment closure method for stochastic reaction networks.
J Chem Phys 2009, 130:134107. PubMed Abstract  Publisher Full Text

Singh A, Hespanha JP: Lognormal moment closures for biochemical reactions.
Proceedings of the 45th Ieee Conference on Decision and Control, Vols 114, IEEE Conference on Decision and Control 2006, 20632068. Publisher Full Text

Gillespie CS: Momentclosure approximations for massaction models.
IET Syst Biol 2009, 3:5258. PubMed Abstract  Publisher Full Text

Konkoli Z: Exact equilibriumstate solution of an intracellular complex formation model: kA ↔ P reaction in a small volume.
Phys Rev E 2010, 82:041922. Publisher Full Text

Mattis DC, Glasser ML: The uses of quantum field theory in diffusionlimited reactions.
Rev Mod Phys 1998, 70:9791001. Publisher Full Text