Kinetic regulation of multi-ligand binding proteins

Background Second messengers, such as calcium, regulate the activity of multisite binding proteins in a concentration-dependent manner. For example, calcium binding has been shown to induce conformational transitions in the calcium-dependent protein calmodulin, under steady state conditions. However, intracellular concentrations of these second messengers are often subject to rapid change. The mechanisms underlying dynamic ligand-dependent regulation of multisite proteins require further elucidation. Results In this study, a computational analysis of multisite protein kinetics in response to rapid changes in ligand concentrations is presented. Two major physiological scenarios are investigated: i) Ligand concentration is abundant and the ligand-multisite protein binding does not affect free ligand concentration, ii) Ligand concentration is of the same order of magnitude as the interacting multisite protein concentration and does not change. Therefore, buffering effects significantly influence the amounts of free ligands. For each of these scenarios the influence of the number of binding sites, the temporal effects on intermediate apo- and fully saturated conformations and the multisite regulatory effects on target proteins are investigated. Conclusions The developed models allow for a novel and accurate interpretation of concentration and pressure jump-dependent kinetic experiments. The presented model makes predictions for the temporal distribution of multisite protein conformations in complex with variable numbers of ligands. Furthermore, it derives the characteristic time and the dynamics for the kinetic responses elicited by a ligand concentration change as a function of ligand concentration and the number of ligand binding sites. Effector proteins regulated by multisite ligand binding are shown to depend on ligand concentration in a highly nonlinear fashion.


Background
A wide variety of intracellular events are initiated via temporal change of ligand concentrations. One of the most important ligands in many cells is calcium (Ca 2+ ). Calcium interacts with and regulates the activities of a large number of calcium-binding proteins as well as numerous effectors. The number of functional calcium binding sites within these proteins can range from one or two to ten or more [1][2][3]. The most common number of calcium binding sites is four: as observed in the most ubiquitous protein, calmodulin as well as troponin and other EF-hand containing proteins [4,5]. Temporal elevation of intracellular free Ca 2+ is the key regulatory factor of the Ca 2+ -dependent protein activity [6][7][8][9][10]. The characteristics of the induced signal are not fully understood. It remains to be determined how a single ligand is able to govern numerous intracellular properties.
Whilst the multisite ligand binding is not limited to Ca 2+ signaling, Ca 2+ is probably the most versatile ion, regulating the largest number of cellular events. Several Ca 2+ -binding proteins can be considered as examples of multisite ligand protein interactions. Structural biology investigations of calcium binding proteins in complexes with target protein peptides have suggested that the specificity in Ca 2+ -CaM binding protein-dependent target activation arises from the diversity of interaction interfaces between the Ca 2+ -regulated protein and its target proteins [5,[11][12][13][14][15][16][17][18][19][20][21]. The most ubiquitous protein, calmodulin (CaM), consists of two globular domains, each domain containing a pair of helix-loop-helix Ca 2 + -binding motifs called EF-hands [1,3,5,16,17]. In earlier studies the authors demonstrated that in addition to the diversity of CaM-target interfaces; the CaM selectivity emerges from its target specific Ca 2+ -affinity; the number of Ca 2+ ions bound and the target specific cooperativity [22][23][24].
Another major factor that contributes to the selectivity of seemingly simultaneous regulation of several multisite Ca 2+ binding proteins and Ca 2+ -mediated processes is the temporal alterations of Ca 2+ [25][26][27]. The remarkable variety of Ca 2+ signals in cells, ranging from infrequent spikes to sustained oscillations and plateaus, requires an understanding of how fast intracellular calcium changes regulate the kinetics of multiple multisite Ca 2+ binding proteins. Therefore mathematical modeling of Ca 2+ jump induced responses could prove to be invaluable in the interpretation of transient kinetic experiments.
Cooperative binding is a special case of molecular interactions where ligand binding to one site of a molecule depends on the ligand binding to the other sites. The first quantitative determination of the dynamic properties of cooperative binding was proposed by [28]. In this work the authors emphasized the significance of the cooperativity by studying the fast dynamics of Ca 2+ binding to calretenin (CR), which has one independent and four cooperative binding sites. The investigation of cooperative effects of Ca 2+ binding to CR was performed both experimentally and using mathematical modeling. The authors employed the simplified version of the Adair-Klotz model [29,30] to describe the dynamics of the interactions involved in Ca 2+ binding to CR. This approach was then extended to the binding of Ca 2+ to CaM [31]. The models proposed in these studies [28,31] demonstrated excellent fitting results to the experimental data, in comparison with the previously published models. However, the described approach is rather limited, as it describes fitting instead of providing a mechanistic description. An alternative methodology offered by [28,31] is not directly applicable from the physical and chemical point of view because the Adair-Klotz model for sequential ligand binding was utilized [29,30] whereas binding of Ca 2+ to EF-hand proteins [32][33][34][35] is non-sequential [22]. Given the importance of studying fast Ca 2+ binding kinetics and the lack of understanding of the underlying mechanisms, we developed a detailed mathematical model for ligand binding to multisite proteins with both cooperative and independent binding sites.
Mathematical modelling of multisite protein kinetics in response to rapid ligand changes presented in this paper provide new insights into the mechanism of conformational kinetics of multisite proteins in complex with variable number of bound ligands for the two distinct physiological situations described: i) When the ligand concentration significantly exceeds protein concentration, ii) When the total amount of ligand is conserved and comparable with the protein concentration.
In the first case, the buffering effects are negligible whereas in the latter, the ligand-protein interactions have a significant impact on the amount of available ligand and the binding kinetics. In this work, the equations for the dependence of the characteristic time constants and the temporal distribution of individual conformations as a function of the ligand concentration, the number of binding sites and the binding affinities have been derived. The impact of the number of binding sites, temporal effects on conformations, and regulation by multisite proteins of their effector proteins have been investigated by employing the developed models. The analysis of the ligand-multisite protein mediated regulation of effector proteins suggests that significant degree of selectivity in regulation can be achieved by a single ligand by employing mechanisms described in this study.

A new model for multisite protein ligand binding kinetics
The majority of studies of the activation of multisite proteins consider only the ligand concentration-dependent profiles. One of the interesting questions about these multisite proteins is how temporal alterations of ligand concentration contribute to their function. The shape of distribution of multisite proteins in complex with variable numbers of bound ligand is known or can be experimentally elucidated in many cases [4,5,16,[22][23][24][36][37][38]. However, the role of temporal transitions caused by fast alteration of ligand concentration on multisite proteins and on multisite protein-regulated target proteins remains unclear. It is reasonable to assume that there can be at least two distinct mechanisms of fast ligand alterationmediated effects exhibited in two distinct system scenarios: i) the ligand concentration is significantly greater than the multisite protein concentration, ii) the ligand is comparable with the multisite protein concentration. In the first scenario ligand binding to multisite protein leads to insignificant changes of free ligand, whereas in the second case, free ligand concentration can vary substantially when ligand molecules bind to the multisite proteins. This paper describes the development of the two models, which address these distinct physiological situations.

The model for abundant ligand concentration
In this model we describe physiological situations where the ligand concentration significantly exceeds the multisite protein concentration. In a previous study [22,24] the authors analyzed functions for the probability of an individual site being in the bound or non-bound state and a function giving the probability of a multisite protein being in a complex with different number of bound ligands (Eqs. (2) and (3) in Methods) [29,39,40]. To investigate the kinetics of the multisite protein ligand interactions the present study extends the previous model to consider the ligand concentration as a function of time (Eqs. (4) and (5) in Methods). The solutions for the individual sites to be in a particular state were obtained for those cases where ligands are subject to rapid changes between steady-states (Eqs. (6) and (7) in Methods). Due to the large number of sites involved, knowledge of the state probability distribution for individual binding sites allows accurate estimation of the dynamics of the total concentration of bound ligand in response to a jump in free ligand concentration (Eq. (11) in Methods).
In order to gain more insights into the distribution of the intermediate protein conformations (complexes with variable number of bound ions) we investigated the case of a multisite protein with identical binding sites (Eq. (12) in Methods). While this case is a relatively rare occurrence in living cells, it enables insight into the role that the number of binding sites plays in cellular signalling. There are several examples of protein families that have variable number of ligand binding sites either due to their structural properties or by them forming large tertiary complexes. For example members of the Ca 2+ family of binding proteins can differ in the number of ligand binding sites [41,42]. The most ubiquitous Ca 2+ -binding protein, calmodulin (CaM), contains four Ca 2+ binding sites as does troponin (TnC) [18] and calcineurin phosphatase (CaN) [43]. However the number of functional Ca 2+ binding sites can vary from two to ten as in the protease, calpain [1] or even more in other cases [2]. To investigate the role that the number of ligand binding sites plays in multisite kinetics, the ligand concentrations, at which the intermediate conformations reach their maximum values, (Eq. (13) in Methods) and the corresponding magnitudes for those conformations (Eq. (14) in Methods) were estimated. Figure 1 shows the maximum protein conformations in complex with one, two and three ligands as a function of the number of binding sites (Eq. (14) in Methods). The graph demonstrates that the magnitude of the ligand-multisite complexes decreases dramatically as the number of binding sites increases. The presented results suggest that the relative magnitude of individual intermediate conformations decreases as the number of binding sites increases. This in turn results in subtler regulatory effects of those proteins with larger number of ligand binding sites. For example, in CaM that has four binding sites for calcium [18], the presence of four sites leads to the increased multifunctionality of this protein due to the additional regulatory properties of intermediate conformations [24]. However, as an exception some forms of CaM have six binding sites [44,45]. In this case, according to Fig. 1 the magnitude of intermediate conformations significantly decreases compared to the case of four binding sites, resulting in decreased regulatory properties of the protein. The presence of more than one binding site results in increased multifunctionality of the protein but at the same time leads the decrease of the regulatory effects. Thus, it seems that there is an "optimal" number of binding sites, which have been developed during the evolution, for instance four calcium binding sites in CaM.
Next, the ligand concentrations for half maximum effective ligand concentrations, U 0 0.5 and U n 0.5 , for the apo-and saturated multisite protein conformations when the protein species equals 50 % of the total concentration were estimated (Eqs. (15) in Methods). This solution shows that the ligand concentration for the half maximal protein activity, known as EC 50 , would be equal to the equilibrium dissociation constant K (EC 50 = K) for proteins with one binding site only (n = 1). Figure 2a shows the dependence of U 0 0.5 /K and U n 0.5 /K, on the number of binding sites. The model predicts that there is a significant change in the required ligand concentration U n 0.5 /K for the fully bound conformation, while U 0 0.5 /K does not change with time. In the previous study the authors reported on the regulatory importance of the distribution of individual multisite protein conformations [22,24]. Here calculations (Eqs. (15) in Methods) are presented for the half-width between the half-maximal effective ligand concentrations U 0.5 /K as a function of the ligand concentration for the intermediate conformations (one bound ligand) of the proteins with different number of binding sites (Fig. 2b) Fig. 2b the half-width range of calcium concentrations is approximately from -1 to 1 on the logarithmic scale, which corresponds to 10 -7 M-10 -5 M due to the fact that the affinity of calcium binding sites in CaM is approximately 10 -6 M [46,47]. Interestingly, the range 10 -7 M-10 -5 M corresponds to the physiological range of intracellular calcium concentrations in cardiac muscle cells [48]. Within this range of calcium concentrations, the switching between calcium channel opening and closure takes place [49].
The difference A n between ligand concentrations U 0.9 /K and U 0.1 /K for the saturated multisite protein conformations, when the protein species are equal to 90 % and 10 % of the total concentration, as a function of the ligand concentration for proteins with different number of binding sites is shown in Fig. 2c and d respectively. The determination of A n allows an understanding as to how an increase of the number of binding sites n affects the steepness of the doseresponse curve (Fig. 2c). It can be seen from Fig. 2c and d that with an increase of n, the width A n decreases while the steepness increases and shifts to the range of higher ligand concentrations U/K. The greater steepness of the dose-response curve caused by the presence of increasing number of binding sites (Fig. 2c) may result in the switch-like response and ultrasensitivity of the protein activation [50]. For instance, the steepness of CaM activation defines the threshold properties for the switching of erythrocyte aggregation and deformability from one steady-state to another [51].
Equation (17) for the total amount of bound ligand in the case of multisite protein with identical binding sites, can be used to estimate the amount of bound ligand when the ligand concentration is equal to the Calculations for the half-width between the half-maximal effective ligand concentrations as a function of the ligand concentration for proteins with two, three, four and five binding sites. c. The difference between ligand concentrations for the saturated multisite protein conformations when the protein species equal to 90 % and 10 % of the total concentration as a function of the ligand concentration for the proteins with one to six binding sites. d. The difference between ligand concentrations for the saturated multisite protein conformations when the protein species equal to 90 % and 10 % of the total concentration as a function of the number of binding sites up to six equilibrium dissociation constant (U = K). Our model predicts that the ligand concentration for the half maximal protein activation, EC 50 , is equal to the equilibrium dissociation constant K for any number of bound sites for a multisite protein with identical binding sites. Figure 3 shows the calculations for the temporal characteristics of the apo-and fully bound species. Figure 3a and b show that the temporal shapes of the apo-and fully bound conformations (Eqs. (19) in Methods) in response to a ligand change are similar to the steady-state dependence of the same conformations on ligand concentration [24]. The kinetic parameters, τ 0 0.5 and τ 4 0.5 can be estimated as the time required to reach 50 % of the total concentration ( Fig. 3a and b). The described kinetic parameters have been investigated as a function of the initial and final ligand concentrations (Fig. 3c, d and Eqs. (24) in Methods). Our analysis reveals a reduction of the time constant, τ 0 0.5 , of the apo-conformation with the reduction of the initial ligand concentration and an increase of the final ligand concentration. However, the dependence of the characteristic time τ 4 0.5 ( Fig. 3d) showed an unexpected bell shaped dependence on the final ligand concentration compared to the simpler monotonic dependence for τ 0 0.5 (Fig. 3c). The model predicts that there is an "optimal" ligand concentration for the saturation effect to take the longest time (Fig. 3d).
The bell shaped dependence of τ 4 0.5 for a protein with four binding sites, for example CaM [18], on final ligand concentration shown in Fig. 3d appears to be related to the presence of intermediate conformations with one, two and three bound sites. In a single site molecule (n = 1), for example crystalline bovine β-trypsin [52] that is characterized by the absence of intermediate forms, τ 4 0.5 depends on U 1 /K monotonically, i.e. there is no bell shaped dependence, as is evident from Eqs.

The affinities of binding sites differentially affect the kinetic responses of intermediate conformations
The proposed model has been employed to investigate the ligand jump-dependent kinetics of both saturated and non-saturated conformations. Initially, an idealized model of a multisite protein with identical binding sites was used to investigate the impact of ligand concentrations on the multisite protein kinetics. However, in living cells there are very few proteins (if any) that have identical ligand binding sites. Therefore the model was extended to examine the implications of variations in binding site affinities on the predicted concentrationresponse profiles. Figure 4 shows the dependence of the time point τ m max k − when the intermediate protein conformations reach the maximum as a function of magnitude of ligand jump (Eqs. (26) in Methods). It can be seen from Eqs.
3 , U 1 /K < 1 and U 1 /K < 3, respectively. Under these special cases, where the ligand concentration U 1 /K is not sufficient for the concentrations of the intermediate conformations to reach their maximal values, these concentrations monotonously grow to their respective steady-state levels. According to Eq. (26), the values U 1 = K < 1 3 ; U 1 /K < 1 and U 1 /K < 3 correspond to the three individual intermediate conformations with one, two and three bound sites respectively.
To investigate the impact of the dissociation constants of individual binding sites we employed the multisite protein model (please see subsection "Multisite proteins with four different ligand binding sites" in Methods) with marginally (Fig. 5) and significantly different association constants (Fig. 6). The main result that follows from the analysis of the intermediate conformation curves is that the affinities of the different binding sites mainly affect the magnitudes of corresponding protein conformation. For example, the conformation of a multisite protein corresponding to the one ligand bound state is present in lower concentration if the affinity of the binding centre is lower. However the overall shape of the concentration dependent profile has not changed. This property is very similar to the case of steady-state dependence on the ligand concentration. The only difference is that the bell shape dependence on time during the kinetic response is partially skewed. However, significant variation in affinities changes the magnitude, and also leads to the asynchronous kinetics of the intermediate conformations.

The effects of cooperativity in Ca 2+ binding to CaM
In order to investigate the influence of cooperativity, we chose a well-characterized protein, CaM, as the model object. The CaM protein contains two independent EFhand globular domains, with two binding sites [1,3,5,16,17]. The sites within each of the domains cooperatively influence each other. It has been reported that cooperative binding occurs between two neighbouring sites within the N-and C-terminal domains of CaM [22,53,54]. Figure 7 shows the model predictions for CaM where we assume that the molecule has two independent domains, with two identical cooperative sites. In the first domain the affinity of one site changes from K 1 = 0.9 μM to K 1 c = 0.2 μM if the other site is occupied and in the second domain the affinity changes from K 2 = 0.8 μM to K 2 c = 0.1 μM (Eq. (28) in Methods) [22]. Figure 7a and b show the influence of cooperativity on the steady-state concentrations of CaM with certain number of bound sites. The presence of cooperativity shifts the dose-response characteristics along the ligand concentration axis and changes the magnitude of intermediate conformations allowing more developed selective regulation of the activity of CaM. The investigation of the dynamic properties of co-operativity in CaM ( Fig. 7c and d) for intermediate, apo-and saturated species revealed that the cooperativity influences the magnitudes of time-dependent characteristics. The proposed model predicts that the cooperative binding leads to more pronounced selective effects for intermediate conformations and higher differences between the initial and steady-state levels for the apo-and saturated forms.
The results of the present analysis suggest that cooperativity plays an important role in the regulation of the activity of multisite proteins by allowing wider possibilities for selectivity. However, the presence of cooperativity leads to  quantitative rather than qualitative changes in the system. The introduction of cooperative binding is crucial for the experimental data fitting but at the same time brings further complexity to the system, which does not necessarily lead to a better understanding of the underlying mechanisms. As a result of this, further analysis was carried out without considering this effect.

The model for comparable ligand and protein concentrations
The previous sections considered the physiological case where the ligand concentration was above saturation level meaning that the ligand-multisite protein interactions did not affect the availability of the ligand. Here the case when the amount of ligand is limited is considered. This can occur in cases when the ligand concentration level is comparable to the multisite protein availability. The model predictions show the relative redistribution of the ligand in the free and bound states. The responses of a multisite protein with four identical binding sites to ligand concentration step change (Eqs. (44) and (45) in Methods) for two different protein concentrations, L T /K = 2 and L T /K = 50, were studied (Fig. 8). Instead of considering absolute ligand concentrations, the approach considered ratios of the ligand concentration to the affinities of the binding sites. A borderline case where the free ligand concentration is barely affected by interaction was considered ( Fig. 8a and b) as well as a smaller total ligand concentration where the free ligand is nearly exhausted as a result of buffering by the multisite protein ( Fig. 8c and d). The comparison of the free ligand concentration (Eqs. (43) and (44) in Methods) for these two cases is shown in   Figure 9a shows the model predictions for the characteristic time required for intermediate protein conformations with one, two and three bound ligands to reach their highest concentrations in response to a step change in ligand concentration (Eq. (50) in Methods). The model predicts that the characteristic time τ 0 0.5 , required for the apo-form to reach its half growth level monotonically decreases with the increase of the total ligand concentration. However, the characteristic time constant τ 4 0.5 , which represents the saturated conformation reveals a distorted bell shaped dependence on ligand concentration ( Fig. 9b and Eqs. (52) in Methods). This bell shaped dependence, which was also observed in Fig. 3d, can also be explained by the presence of intermediate conformations. The distortion of the bell shape in Fig. 9b appears to be due to the ligand consumption that is included into the consideration in this section and was not considered in Fig. 3d. This result may be significant for the dynamics of CaM activation as our model predicts that with an increase of the total ligand concentration, the limited amount of ligand leads to an additional increase of τ 4 0.5 compared to the case without the ligand consumption. The model, therefore, predicts possible transient differences in multisite protein signal transduction in response to fast transient kinetics of multisite proteins. ). The model predicts that due to the lack of available ligand and buffering by the multisite protein in the case of limited amount of ligand, the multisite protein is unable to become fully saturated after the step change in ligand, and the majority of the ligand becomes distributed among the intermediate species. e. The comparison of the dynamics of free ligand concentration U/K after step change in ligand. The amount of available ligand is barely altered for L T /K = 2, and exhausted when the ratio of total protein concentration to the binding constant is L T /K = 50

Discussion
In this paper we analysed multisite protein kinetics in response to rapid changes in ligand concentrations. The model for multisite protein kinetics for variable number of bound ligands was developed for two physiological cases: when the concentration of ligand is much higher [55,56] and when it is comparable with the concentration of the multisite protein [57]. The results obtained by the proposed model allow for an accurate interpretation of the experimental data for the concentration of multisite proteins such as CaM, TnC, CaN and other Ca 2+ dependent secondary messengers regulated by Ca 2+ ions. The kinetic effects in response to ligand binding [13,58,59], required for understanding of pathway regulation, can be interpreted by the presented model.
The approach developed in this project is applicable to a number of areas of kinetic experiments. One of them is widely used techniques to study chemical kinetics using the pressure jump technique [12,60,61]. According to our model, the effects induced by rapid change in pressure leading to the change in protein-ligand interactions [12,60,61], can be interpreted and explained by the alteration in the affinities of the binding constants of multisite proteins. The time dynamics of the individual multisite protein species can offer new insights into the underlying biophysical mechanism of ligand-protein interactions in response to fast change.
Our model shows that the concentration of intermediate conformations as a function of time represents skewed bell shapes. We found that as the protein concentration rises, free ligand concentration becomes exhausted [13]. This result is consistent with experimentally observed Ca 2+ -CaM-dependent inhibition of myosin motor function [62]. The results obtained by this model increase our understanding of differential activation of protein phosphatase 2B (PP2B) [59,63] and calcium/calmodulin-dependent protein kinase II (CaMKII) [64,65] kinetics. PP2B binding increases the affinity of CaM for its targets [59] and, therefore, is likely be activated by low amounts of calcium.
The presented model offers a new tool for the interpretation of transient kinetics experiments performed by the flash photolysis and stopped-flow techniques [66][67][68]. Availability of a caged analogue of the necessary reactant is considered as one of the biggest problems of flash photolysis [66], which can be overcome by employing the proposed model with the low amount of protein. A less obvious area of application for this methodology is the kinetics of proteins with multiple phosphorylation sites [69][70][71][72]. This work shows that the highly versatile intracellular multifunctionality of multisite proteins is achieved not only by the order of ligand-protein interaction and the number of bound ligands, but also by temporal regulation.

Conclusions
The developed in this study models for the kinetics of the multisite ligand-receptor binding for the two physiological cases, where the ligand concentration is abundant and comparable with the protein concentration, make universal predictions for the temporal distribution of multisite protein conformations in complex with variable numbers of ligands. The strongest effect of the ligand availability is observed for the multisite protein conformations with larger numbers of bound ligands. The two models show that the concentration of individual multisite protein conformations changes with time nonlinearly and that the temporal distribution for the concentrations of the intermediate conformations represents skewed bell shapes. The models derive the characteristic times and the dynamics for the kinetic responses elicited by a ligand concentration change as a function of ligand concentration and the number of ligand binding sites. The developed models allow for a novel and accurate interpretation of concentration and pressure jumpdependent kinetic experiments. Our models are applied to study the kinetics of calmodulin, however the models also provide universal predictions and allow us to extend the understanding of a large number of multisite bindingregulated circuits and the mechanisms underlying dynamic ligand-dependent regulation of multisite proteins.

Kinetic properties of multisite proteins with n binding sites
Here, we describe mathematical equations used to describe dynamic interactions of ligand molecules with multisite proteins. The model described in this section extends our previous analysis for multisite protein interactions under steady-state conditions [22][23][24].

Multisite protein with independent ligand binding sites
The kinetic scheme for such interaction can be represented as follows: where L i 0 is the ith binding site of the multisite protein in unbound state, U is a ligand molecule, L i 1 is the ith binding site of the multisite protein being occupied, k i + and k i − are the association and dissociation rates respectively.
The probabilities for the ith binding site to be not occupied or occupied as a function of ligand concentration U are given by: where There are two possible states of a binding site: occupied or not occupied. Since the number of sites in the molecule is n, there are 2 n possible molecular forms, i.e. the states characterized by combinations of bound and free sites.
The probability for a multisite protein with independent binding sites to be in a particular molecular form is given by multiplications of probabilities of ligand binding at each site: where j ¼

The kinetics of multisite protein interactions with abundant ligand concentrations (U T > > L T )
In this case we assume the ligand concentration only changes at a given time point (t = 0) in a step fashion from U 0 to U 1 and remains constant otherwise.
The probability of a multisite protein being occupied by ligand molecules as it was shown in Eq. (3) as a function of time can be reformulated as follows: Þ can be determined by considering the original set of ordinary differential equations for single site interaction of a protein with a ligand according to Eq. (1) [24]: where L i 0 (U, t) is the concentration of a free binding site, L i 1 (U, t) is the concentration of a bound ligand molecule and L T is the total number of the protein molecules.
Here we use steady-state solutions K i þU as initial conditions for the ligand concentration jump from U 0 to U 1 . A particular solution for the system of differential Eqs. (5) in response to the ligand concentration shift from U 0 to U 1 is given by: where τ U 1 ð Þ ¼ The normalisation of the solution (6) by the total protein concentration allows the definition of probability of the ith binding site to be in an occupied p i 1 (U, t) or unoccupied p i 0 (U, t) state, respectively, at a given ligand concentration: At any given time the total concentration of the protein, L T , is conserved and the sum of the probabilities (7) equals 1: In the most general case, the concentration of the jth molecular form, M j (U, t), of a molecule with n different binding sites in response to a step in ligand concentration is given by the product of probabilities (7) according to Eq. (4): where c i (j) = 0 or c i (j) = 1 for free or occupied ith binding site respectively. The probabilities for individual molecular forms in steady-state can be obtained from Eq. (9) by setting t → ∞: The kinetics of the amount of ligand bound to a multisite protein with n binding sites can be written as follows: Multisite proteins with identical ligand binding sites The multisite proteins that contain n identical binding sites with equilibrium dissociation constant equal K ¼ k − k þ can then be considered. The steady-state concentration of one of the molecular forms of multisite protein with n identical binding sites for m bound sites, as a function of ligand concentration [22,24] is given by: For intermediate forms the N 1 (U), …, N n − 1 (U) species dependence on ligand concentration represents a bell shape with one molecular form "magnitude" at a particular ligand concentration U m max . Differentiating Eq. (12) with respect to U and solving dN/dU = 0 for U gives: The magnitudes N m max of the intermediate conformations N m corresponding to U m max values are: The multisite protein conformations in the apo-, N 0 , and in the fully saturated states, N n , would reach their maximum that equal to the total multisite protein concentration L T under conditions of very low and very high ligand concentrations, respectively.
Equation (12) can be used to estimate the half maximal effective ligand concentration (EC 50 ), U 0 0.5 and U n 0.5 , for the apo-and saturated multisite protein conformations respectively, when the protein species equal 50 % of the total concentration L T : It can be noted that the half maximal effective ligand concentration is equal to equilibrium dissociation constant (EC 50 = K) in proteins with only one binding site (n = 1).
Equations (15) can be solved with respect to n: Equation (11) can be used to derive the full amount of ligand bound to multisite protein for U = K: where U sat is the ligand concentration for the case when all binding sites are occupied. The amount of bound ligand for U = U sat is given by S(U sat , t) = n ⋅ L T and for U = K is given by The dynamic alterations of intermediate conformation N m (U, t) in response to ligand concentration jump from U 0 to U 1 according to Eqs. (6) and (7) are given by: Differentiating Eq. (19) with respect to t and solving for dN m /dt = 0 yields the time τ m max when the multisite protein forms N m with m bound molecules of ligand reach their maximal values N m max : The substitution of τ m max into Eq. (19) gives the maximal values of intermediate multisite protein conformations reached at τ m max : Comparison of Eqs. (14) and (21) suggests that the steady-state maximum values of intermediate multisite protein conformations steady-state are similar to those transiently reached during the dynamic response to the step in ligand concentration.
According to Eq. (19) N m (U, t) for the apo-and fully saturated forms when t = 0: According to Eq. (12) steady state levels for the apoand fully saturated forms when t → ∞: Equation (19) is further used to define the time,

Multisite proteins with four identical ligand binding sites
This study next considered the kinetic properties of a protein with 4 binding sites. This allows 2 4 = 16 molecular forms, each with potentially unique biochemical properties. There are 4 possible combinations of protein species bound to one or to three ligands, and 6 possible distinct molecular forms with two sites occupied. In the previous work we described the steady-state dependence of the individual multisite conformations on ligand concentration [22][23][24]. Here we analyse the kinetic transition of the individual species concentrations in response to the step in ligand concentration.
The dynamical alterations of intermediate conformation N m (U, t) in response to a step in ligand concentration from U 0 to U 1 according to Eq. (19) are given by: where: are the dissociation and equilibrium dissociation constants for ligand binding, respectively.
The maximum values of N m max are reached at the following ligand concentrations: u max , for the multisite protein species with one, two and three bound ions, respectively according to Eq. (13) and equal N 1 max = 0.105L T , N 2 max = 0.063L T , N 3 max = 0.105L T , according to Eq. (14).
Differentiating Eq. (25) with respect to η and solving dN m /dη = 0 for η gives the non-dimensional time η m max = τ m max k − when the intermediate species reach their maximum: Equations (26) suggest that η 1 max , η 2 max and η 3 max are undefined when u 1 < 1 3 ; u 1 < 1 and u 1 < 3; respectively (q 1 , q 2 and q 3 asymptotes). Under these conditions, the intermediate species do not reach the maximum in response to ligand step, instead their relative number increase in a monotonous manner.

Multisite proteins with four different ligand binding sites
In this section we analyse the kinetic properties of a multisite protein with four different binding sites. The steady-state analysis can be found in the authors previous investigation [22,24].
It is assumed that all association k 1 + , k 2 + , k 3 + , k 4 + and dissociation , k 4 − rates are unique for each binding centre. Then assuming for example that k 1 The non-dimensional concentration u ¼ U K and non-dimensional constants h 1 = 1, The set of Eqs. (25) can then be employed to calculate the dynamics of the multisite protein conformations bound to a different number of ligand molecules.

Multisite proteins with two pairs of cooperative binding sites
The molecule contains two independent domains A and B, with two identical cooperative binding sites. The domain A is described as follows: The ODEs for the scheme (27) is given by: The total number of species that follows from Eqs. (28) is given by: The steady-state solutions of the system (28) are given by: where One can re-write Eqs. (30) as follows: are the probabilities for the domain to be in a particular conformation due to the bound ligand molecules.
The probabilities for the other domain, which also contains a pair of cooperative binding sites, are given by: These probabilities were derived for the two domains of molecule (Eqs. (31) and (32) respectively). The probabilities of the molecule to be in a certain conformation with 0, 1 or 2 bound ligands in each of the domains, are as follows: where p i,j (U) is the probability of the protein conformation with i bound sites in the first domain and j bound sites in the other. We use the sum of probabilities for the case of 1 bound site in a domain since we assume that all the sites in each domain are identical. The concentrations of the molecular forms of the protein with certain number of bound sites are given by: One can rewrite Eq. (34) as follows: Next we find the kinetic solution of system (28) for U = U 1 : The kinetics of multisite protein interactions with constant ligand concentrations In this section we consider the mechanism of multisite binding for the case when the total ligand concentration is conserved. Under this assumption the system of differential equations for ligand binding to a molecule with single binding site (5) needs to be complemented by the law of ligand conservation: where L 0 is the concentration of the free site, L 1 is the concentration of the occupied site, U T and L T are the total concentrations of ligand and protein molecules, respectively. There is only one positive steady-state solution of the system (37): where K ¼ k − k þ and F U T We use the steady-state solutions (38) as initial conditions to find the particular solution. The solution of the system (37) in response to the ligand concentration jump from U T0 to U T1 is given by: where C U T ð Þ ¼ The normalisation of the solution (39) by the total protein concentration allows the definition of probability of the binding site to be in an occupied p 1 (U T , t) or unoccupied p 0 (U T , t) state, respectively, at a given total ligand concentration: where η = t ⋅ k − . The concentration of the jth molecular form, M j (U T , t), of a protein with n independent binding sites in response to a step change in ligand concentration is given by the product of probabilities according to Eq. (9): where c i (j) equals 0 or 1 for free and occupied sites, respectively and j ¼ The kinetics of the amount of ligand bound to a multisite protein with n independent binding sites can be written as follows: The concentration of free ligand can be written as the difference between the total ligand concentration and the bound ligand concentration (42): The probabilities of the ith site to be free or occupied respectively for the molecule with n independent binding sites in this case: where F U T ð Þ ¼ The dynamic alterations of the molecular forms with m bound out of n independent identical binding sites N m (U T , t) in response to the total ligand concentration jump from U T0 to U T1 according to Eqs. (19) and (44) are given by: The concentration of multisite protein conformations, N m , bound to m ligand molecules as a function of ligand concentration in steady-state is given by: The ligand concentrations, U m max , for the maximal values of intermediate protein conformations, N m (U m max ), can be found by differentiating Eq. (46) with respect to U T and solving dN m /dU T = 0 for U T :