Resolving the Fast Kinetics of Cooperative Binding: Ca2+ Buffering by Calretinin

Cooperativity is one of the most important properties of molecular interactions in biological systems. It is the ability to influence ligand binding at one site of a macromolecule by previous ligand binding at another site of the same molecule. As a consequence, the affinity of the macromolecule for the ligand is either decreased (negative cooperativity) or increased (positive cooperativity). Over the last 100 years, O2 binding to hemoglobin has served as the paradigm for cooperative ligand binding and allosteric modulation, and four practical models were developed to quantitatively describe the mechanism: the Hill, the Adair-Klotz, the Monod-Wyman-Changeux, and the Koshland-Némethy-Filmer models. The predictions of these models apply under static conditions when the binding reactions are at equilibrium. However, in a physiological setting, e.g., inside a cell, the timing and dynamics of the binding events are essential. Hence, it is necessary to determine the dynamic properties of cooperative binding to fully understand the physiological implications of cooperativity. To date, the Monod-Wyman-Changeux model was applied to determine the kinetics of cooperative binding to biologically active molecules. In this model, cooperativity is established by postulating two allosteric isoforms with different binding properties. However, these studies were limited to special cases, where transition rates between allosteric isoforms are much slower than the binding rates or where binding and unbinding rates could be measured independently. For all other cases, the complex mathematical description precludes straightforward interpretations. Here, we report on calculating for the first time the fast dynamics of a cooperative binding process, the binding of Ca2+ to calretinin. Calretinin is a Ca2+-binding protein with four cooperative binding sites and one independent binding site. The Ca2+ binding to calretinin was assessed by measuring the decay of free Ca2+ using a fast fluorescent Ca2+ indicator following rapid (<50-μs rise time) Ca2+ concentration jumps induced by uncaging Ca2+ from DM-nitrophen. To unravel the kinetics of cooperative binding, we devised several approaches based on known cooperative binding models, resulting in a novel and relatively simple model. This model revealed unexpected and highly specific nonlinear properties of cellular Ca2+ regulation by calretinin. The association rate of Ca2+ with calretinin speeds up as the free Ca2+ concentration increases from cytoplasmic resting conditions (∼100 nM) to approximately 1 μM. As a consequence, the Ca2+ buffering speed of calretinin highly depends on the prevailing Ca2+ concentration prior to a perturbation. In addition to providing a novel mode of action of cellular Ca2+ buffering, our model extends the analysis of cooperativity beyond the static steady-state condition, providing a powerful tool for the investigation of the dynamics and functional significance of cooperative binding in general.


Introduction
In all eukaryotic cells, Ca 2þ signals play a crucial role in the regulation of many cellular processes, including gene expression, cytoskeleton dynamics, cell cycle, cell death, neurotransmission, and signal transduction. To achieve its role as messenger, the intracellular Ca 2þ concentration ([Ca 2þ ]) is very tightly regulated in time, space, and magnitude. The spatiotemporal characteristics of short-lived and often highly localized changes in intracellular [Ca 2þ ] result from a complex interplay between Ca 2þ influx/ extrusion systems, mobile/stationary Ca 2þ -binding proteins (CaBPs), and intracellular sequestering mechanisms. Understanding the kinetics of cellular Ca 2þ transients and its influence on Ca 2þ -regulated processes requires a precise knowledge of the Ca 2þ sensitivities and binding properties of all the components involved, including the binding dynamics to buffering and signaling CaBPs. However, uncertainties in current models studying intracellular Ca 2þ signaling arise mostly from the lack of accurate data on the binding properties of specific molecules involved in Ca 2þ handling, considerably limiting the value of such modeling [1]. An important step towards the goal of precisely describing intracellular Ca 2þ transients was the study by Nagerl et al. [2], in which the relevant parameters (affinities and on-and offrates of Ca 2þ binding) for the CaBP calbindin D-28k (CB) were determined in vitro by flash photolysis of caged Ca 2þ . Cooperative binding of Ca 2þ , known to play a significant role in multisite CaBPs such as calmodulin [3] and calretinin (CR) [4], has never been directly determined in rapid kinetic experiments, but only inferred from steady-state conditions using Hill [5] and Adair-Klotz models [6,7].
Cooperativity first evidenced by oxygen binding to hemoglobin [8] is considered one of the most imperative functional properties of molecular interactions in biological systems, even considered to be the great secret of life, second only to the structure of DNA [9]. Cooperativity is the ability to influence ligand binding at a site of a macromolecule by previous ligand binding to another site of the same macromolecule. Many proteins show increased (positive cooperativity) or decreased (negative cooperativity) affinity for a ligand after binding of a first ligand. Over the last 100 years, hemoglobin has been a paradigm for cooperative ligand binding and allostery. Oxygen binding to hemoglobin resulted in four commonly used descriptions for cooperativity (for review see [10]): the Hill [11], the Adair-Klotz [6,7], the Monod-Wyman-Changeux (MWC) [12], and the Koshlan-Némethy-Filmer (KNF) [13] models. Yet all these models describe cooperativity only when the binding reactions are at equilibrium. Since temporal aspects of most ligand binding processes are essential for correct physiological functioning, it is imperative to consider the kinetics of cooperative binding. To date, studies determining the kinetics of cooperative binding to biologically active molecules have been carried out using the MWC model, in which cooperativity is established by assuming two allosteric isoforms with different binding properties. These studies were limited to special cases where transition rates between allosteric isoforms are much slower than the binding rates [14,15] or where binding and unbinding rates could be measured independently [16]. For all other cases, the mathematical description becomes too complex for simple interpretations [10,17]. The Hill equation is perhaps the oldest and most widely used description for the relative amount of binding by a cooperative molecule, and the cooperative binding is described with two constants: the dissociation constant (K d ) reports on the concentration of ligand at which the cooperative molecule is half occupied and the Hill number (n H ) describes the steepness of the binding curve at the value of K d , denoting a simple quantification of cooperativity. Although not representing a mathematically correct description of cooperative binding at equilibrium-a fact that is stated in the original work [11]-the Hill equation has proven to be extremely useful, as it describes occupancy as a function of ligand concentration with merely two constants that are easy to interpret intuitively. With this in mind, we wanted to resolve the kinetics of Ca 2þ binding to CR and to find an intuitively ''accessible'' quantitative description of the binding kinetics.
CR belongs to the superfamily of EF-hand Ca 2þ -binding proteins. This superfamily is named after the common Ca 2þ binding structure-the EF-hand-first described as the Cterminal E-helix-loop-F-helix Ca 2þ binding site in parvalbumin [18]. Most members have an even number of EF-hand domains organized in pairs [19], representing a structurally conserved architectural unit. CR has six EF-hand domains [4,[20][21][22][23][24], which can be subdivided into two independent domains: one with the cooperative pair of binding sites I and II, and another with binding sites III-VI [23]. Sites III-VI can be further subdivided into one cooperative pair, sites III and IV and sites V and VI [4]. Of the latter pair, only site V binds Ca 2þ , whereas site VI is ''inactive'' [4]. Thus, CR has two pairs of cooperative binding sites (I-II and III-IV) and one independent binding site (V) for Ca 2þ . We tried several approaches to describe the kinetics of these five binding sites based on the published models, and discovered a new and simplified kinetic model that quantitatively resolves the kinetics of cooperative binding. This new model also revealed unexpected and highly specific nonlinear properties of cellular Ca 2þ regulation by CR.

The Binding Sites of Calretinin
We determined the kinetics of Ca 2þ binding to CR by Ca 2þ uncaging using DM-nitrophen (DMn) and measuring the changes in [Ca 2þ ] with the fluorescent Ca 2þ indicator dye Oregon Green BAPTA-5N (OGB-5N) as previously described [2,25]. Changes in the OGB-5N fluorescence were observed immediately after photolysis of DMn in solutions containing various concentrations of CR ( Figure 1A). A rapid rise in [Ca 2þ ] ensued from a resting concentration of approximately 2.4 lM to 11 lM. At this initial [Ca 2þ ], approximately 99.5% of the DMn is in the Ca 2þ -bound form, ensuring that (1) virtually every uncaged DMn molecule will release Ca 2þ , and (2) the amount of free DMn capable of rebinding uncaged Ca 2þ ions [25,26] is considerably limited. This is evidenced by the negligible drop in [Ca 2þ ] in the absence of CR, leading to an almost step-wise increase in [Ca 2þ ] ( Figure 1A). The presence of CR (31 and 62 lM) resulted in a [CR]-dependent drop in [Ca 2þ ]. Figure 1B depicts a simplified reaction scheme of the experiment. To determine the association and dissociation rates of Ca 2þ binding to CR, all these different reactions were incorporated into a mathematical model (see below). The aim was to find a mathematical description for the Ca 2þ -binding properties of CR that best fits the experimental fluorescence traces generated under various conditions. As a starting point to determine CR's kinetics of Ca 2þ binding, we relied on steady-state Ca 2þ -binding properties determined previously. With the same human recombinant CR, selected Hummel-Dryer experiments yielded a Hill

Author Summary
The binding of a ligand to a protein is one of the most important steps in determining the function of these two interactive biological partners. In many cases, successive binding steps occur at multiple sites such that binding at one site influences ligand binding at other sites. This concept is called cooperative binding, and constitutes one of the most fundamental properties of biological interactions. The functional consequences of cooperativity can be accurately resolved when reactions are at equilibrium, but mathematical complexity has prevented insights into the dynamics of the interactions. We studied the protein calretinin, which binds Ca 2þ in a cooperative manner and plays an important role in shaping Ca 2þ signals in various cells. We used two models, a widely tested one and a novel, mathematically simplified one, to resolve the dynamics of a cooperative binding process. The cooperative nature of Ca 2þ binding to calretinin results in accelerated binding as calretinin binds more Ca 2þ . This behavior constitutes an important new insight into the regulation of intracellular Ca 2þ that cannot be matched by noncooperative artificial Ca 2þ buffers. Our simple mathematical model can be used as a tool in determining the kinetics of other biologically important molecular interactions. coefficient of 1.3 for the four binding sites, with a K d of 1.5 lM [4]. However, by flow dialysis, the steady-state binding of Ca 2þ to human CR could be described with the following macroscopic constants (K 1 through K 5 ) 2.2 3 10 5 M À1 , 3.2 3 10 5 M À1 , 4.7 3 10 5 M À1 , 8.0 3 10 5 M À1 , and 2.0 3 10 4 M À1 [4]. The resulting binding curve derived from these values could be accurately fitted with two Hill equations; one equation described four cooperative binding sites with a K d of 2.5 lM and a Hill coefficient of 2.4, and the other one described a single independent site with a K d of 53 lM. In agreement with this data, equilibrium dialysis experiments with chick CR revealed a Hill coefficient of 1.9 [22]. Even Hill coefficient values of up to 3.7 for Ca 2þ -induced tryptophan (Trp) fluorescence changes in rat CR have been reported [21]. However, these conformational changes measured by Trp fluorescence do probably not linearly relate to Ca 2þ binding. Thus, the absolute values should be interpreted with caution, but nonetheless, cooperativity of Ca 2þ binding also occurs in rat CR. Based on these Hill coefficients of 1.3, 1.9, and 2.4, we conjectured that CR has four binding sites for Ca 2þ with positive cooperativity, with a Hill coefficient of approximately 2 and one independent binding site for Ca 2þ .
In accordance with these and other earlier findings on the structure and physiology of CR (see Introduction), we modeled the protein as possessing two pairs of cooperative binding sites (B I B II ) and (B III B IV ), and one independent binding site B V . Because the properties of the two cooperative pairs in CR were considered indistinguishable in the steady-state study [4], thus indicating that the cooperative binding sites are fairly similar, we assumed the properties of both cooperative pairs to be identical: Such an assumption is most useful for reducing the number of variables, thus increasing the reliability of fitting procedures by constraining the model.

Analysis of Ca 2þ Binding to CR Using a New Model That Includes Cooperativity
The most straightforward approach to determine the kinetics of a system would consist of fitting the [Ca 2þ ] decay with a set of exponential functions. However, rebinding of Ca 2þ to free DMn affects the decay kinetics. In addition, the changing properties of the binding sites that underlie cooperativity are expected to cause a shift in the kinetic properties of binding during the decay phase. As a consequence, the relative contribution of multiple decay time constants is continuously shifting. Although fitting the [Ca 2þ ] decay with exponentials might result in time constants for a given trace, this does not allow accurate deduction of the Ca 2þ -binding kinetics of CR. A mathematical model simultaneously describing all processes taking place in the recording chamber, including a ''total'' description of the cooperative and noncooperative binding of Ca 2þ to CR is expected to yield more reliable information on the kinetic properties of CR. To model the cooperativity, we started out by including an allosteric influence between the binding sites of the pairs. This was achieved by setting two states (R and T) for a particular binding site, with each having its own set of rate constants. A binding site is in the ''tensed'' state (T), with a low affinity for Ca 2þ , when no Ca 2þ is bound to the other site in the pair, whereas a binding site is in the ''relaxed'' state (R), with a high affinity for Ca 2þ , when the other site already has a Ca 2þ ion bound. We assumed that binding of Ca 2þ to one site always leads to a rapid transition T!R in the other site and that an unbinding of Ca 2þ from one site always leads to a rapid transition R!T in the other site. This allowed us to incorporate the transition rates between states R and T in the binding and unbinding rate constants, further simplifying the model (Figure 2A). CR was thus modeled as if consisting of two independent proteins described by the following binding reactions: This cooperative part of the model can be easily related (see Discussion) to the Adair model [6], which provides the most general description of equilibria in terms of stochio- metric binding. For the independent site of CR, we used a standard equilibrium equation: where k on(R or T) and k off(R or T) are the association and dissociation rate constants for the individual cooperative binding sites depending on their Ca 2þ -binding status, and k on(V) and k off(V) are the rate constant for the independent site. The total concentrations of the different ''virtual'' parts are: Despite the simplifying assumptions concerning the cooperative sites, the model that allows a fitting routine to proceed is fairly complex, because it has a considerable number of degrees of freedom. Thus, a procedure was developed that significantly constrains the fit to minimize the variance of the fit results. Simultaneously fitting combined sets of uncaging data (Figure 3; see Materials and Methods) obtained under different experimental conditions sufficiently constrained the model to yield consistent results. We performed a number of individual uncaging experiments generated at one of seven initial conditions A-G (Table in Figure 3). These conditions varied in the initial free [Ca 2þ ] (hence, total [Ca 2þ ]), total [CR], total [DMn], as well as on the lot number of OGB-5N, with each lot having slightly different properties (see Materials and Methods). Under each con-dition (A-G), we performed 12 to 25 uncaging experiments, each one with a different flash energy of the UV laser leading to different amounts of uncaged Ca 2þ . In total, 123 traces were obtained, covering a wide rage of uncaging energies and, subsequently, a wide range of increases in [Ca 2þ ] (see gray areas in Figure 3; for all 123 individual traces, see Figure S1). We set out to find a satisfactory model that would be able to fit all curves obtained under conditions A-G and all tested uncaging intensities. The obtained results from the modeling should be able to describe all experimental curves with a unique set of parameters describing the kinetic properties of CR. To confine the fits, 38 sets of 14 pseudo-randomly picked traces consisting of two traces from each initial condition A-G were generated. The 38 sets were chosen randomly, with the precondition that every trace of a specific starting condition A-G was represented equally. To create 38 sets, each individual measurement was picked at least three and at most eight times. On average, each trace was picked 4.3 times (for details, see Figure S2). Each set of traces was fitted with the model, and the fitted parameters describing the properties of CR were constrained to be identical for all individual traces within one set. The only variable parameter between traces was the amount of uncaging that was fitted individually for each trace. The new model was programmed as to fit the K d(V) and k on(V) for the independent site. However, to aid the choice of starting values, the cooperative part of the model was set up such that k on(R) , k on(T) , the apparent K d (K d(app) ) for the pairs, and the Hill number (n H ) could be fitted. This was achieved by adding a calculation step that determined K d(R) and K d(T) from the latter two parameters (see Protocol S1): Previously determined steady-state parameters (apparent) K d 's and n H of CR [4,[20][21][22] served as starting points in the modeling and helped to further constrain the model. Various combinations of k on starting values between 10 5 and 10 8 M À1 s À1 were tested, but this did not significantly influence the outcome of the fit, indicative of the ''robustness'' of the modeling procedure. Occasionally a particular set out of the 38 yielded an atypical fit with values significantly deviating from the general population of results. If this was the case, we followed up with two approaches. First, we tested whether any of the other 37 sets of fluorescence traces could also be fitted with these deviating values, which in almost all cases yielded unsatisfactory fits. Second, we tested whether the deviating set of traces could also be fitted with the more homogeneous values of the general population of sets by choosing starting parameters closer to these values. In this case, the deviating set could always be fitted with values comparable to the homogeneous constants. It should be noted that the critical parameters (K d(V), k on(V), k on(R) , k on(T) , K d(app) , and n H ) were never constrained or fixed to a certain value. The atypical fit Figure 2. Models Used to Simulate Cooperativity Both models used to simulate cooperativity are based on the ability of the binding sites to occur in two states, one with a low affinity for Ca 2þ (T) and one with a high affinity (R). In the new model (A), a binding site is in the T state when the other binding site has no Ca 2þ bound, and it is in the R state when the other site has Ca 2þ bound. In this model, the transitions from T to R and vice versa are considered to occur instantaneously. In the MWC model, the cooperative pair of binding sites are either both in the T state or both in the R state. The transitions between R and T occur according to the transition constants k þ and k À . For each occupation level, there is a separate equilibrium (L 0 , L 1 , and L 2 ) for which it can be shown that under specific conditions, the equilibrium shifts from mainly the T state in the unoccupied cooperative sets to mainly the R state in fully occupied states (see also Equations 8,9,and 10). doi:10.1371/journal.pbio.0050311.g002 results were probably caused by local minima in the error function of the fit routine. Such local minima are expected, based on the fairly large number of degrees of freedom where the parameters are not completely independent. Deviations in one parameter can be partially compensated by ''shifting'' other parameters. Initially, we used a model that did not include cooperativity and found that most of the 38 individual sets could be fit reasonably well. For a given individual dataset, the quality of the fits were similar between a model with a Hill coefficient of either 1 or 1.9. But when comparing the fit results of all 38 sets, most of the fitted parameters showed strong deviations (up to five orders of magnitude, depending on how the model was exactly defined) when using n H ¼ 1. The high variability of the binding parameters found when assuming n H ¼ 1 is shown in Figure  S3. This indicated that there is no unique solution to describe CR's Ca 2þ -binding properties without cooperativity, in line with previous steady-state findings of n H values between 1.3 and 2.4 [4,22]. Thus, only when including cooperativity and starting the modeling procedure with previously determined steady-state parameters for CR [4,[20][21][22]] did we find a congruent set of values for the fitted parameters for all 38 sets of 14 traces.
An accurate model describing CR's Ca 2þ -binding dynamics should be able to fit all the experimental traces obtained under any condition. This should be the case at lower resting [Ca 2þ ], when Ca 2þ -free binding sites of any affinity in any state are abundant ( Figure 3A-3D and 3G), but also at higher [Ca 2þ ], when mostly the lower affinity independent site V is available ( Figure 3E and 3F). Furthermore, the model is also able to closely describe the [Ca 2þ ] signals after a relatively large uncaging, when [Ca 2þ ] is so high that the buffering by CR is relatively small (see upper trace in Figure 3F). The goodness of fit of the fit procedure can be appreciated by the averaged error for the 38 fits ( Figure 3H, black bars), which shows systematic deviations (if any). Errors were found to be extremely small for the first 20 ms after the flash, and at time points greater than 20 ms, the averaged fits show a small systematic undershoot of the experimental data, yet never exceeding 1.5% of the actual amplitude ( Figure 3H, black bars). The larger errors towards the end of the traces are likely due to the small amplitudes of the signals at these time points, which increase the relative error when there is a constant absolute deviation. But even the largest deviations of the 38 fits (the striped bars showing the largest deviation in either direction) never exceeded 5% of the measured amplitude. The average absolute error (not allowing for positive and negative errors to cancel each other) was maximally 2.1% and again only found towards the end of the traces. Thus, the fit procedure, applying our model, allowed accurate quantifying of the kinetic properties of CR.
All 38 results for the fitted values were plotted as a lognormal cumulative probability distribution because they have a log-normal distribution (Figures 4 and 5), except for the n H value, which was normally distributed ( Figure 5E). These results were then fitted with a normal distribution to determine average and standard deviation. The results of these fits are shown in Table 1. To summarize, we conclude that CR can be described with one independent binding site with a K d of 36 lM, a k on of 7.3310 6 M À1 s À1 , and a k off of 240 s À1 together with two identical cooperative pairs of binding sites with an initial (T state) K d of 28 lM with a k on of 1.8310 6 M À1 s À1 and a k off of 53 s À1 that will dramatically change to a (R state) K d of 68 nM with a k on of 3.1 3 10 8 M À1 s À1 and a k off of 20 s À1 once the cooperative partner site has already bound Ca 2þ . The K d(app) for the cooperative sites is 1.4 lM with an n H of 1.9.

Modeling Cooperativity of Ca 2þ Binding to CR Using a Monod-Wyman-Changeux Model
To compare the results obtained with our new ''simple'' model, we subjected our experimental data to a previously described cooperative model, the MWC model. There, the protein as a whole can switch between two states: one state in which all the binding sites are in the T state, and another one in which all the binding sites are in the R state ( Figure 2B).
The equilibrium between the T and R states is described with the equilibrium constant L, which is also dependent on the number of Ca 2þ ions bound (see Figure 2B). It can be shown [12,27] that: where k þ and k À are rates of the transition between the two states R and T. The indices (0 and i in Equations 8 and 9, respectively) indicate the number of Ca 2þ ions bound. K is often indiscriminately used for both association and dissociation equilibrium constants; here, we denote K as association equilibrium constants, whereas we use K d for dissociation equilibrium constants. Furthermore, the pair of cooperative binding sites can transition from R to T and back, independently of the number of sites that are occupied, thus all transitions defined by L 0 , L 1 , and L 2 are possible. However, for steady-state purposes, one transition (L 0 ) suffices, as described in the original MWC model [12]. Although steadystate properties are independent of the transitions allowed, the kinetic properties will highly depend on the number of allowed transitions. We chose to allow all possible transitions because it was used in earlier kinetic fits with this model [14,16]. We also attempted to fit the data with a MWC model in which only the L 0 transition was possible; however, the resulting fits showed large deviations (.20%) and generated traces with significant deviations from the experimental data.
From Equation 9, we can derive the information that while K T , , K R and L 0 . . 1, the equilibrium between the T and R states is shifted towards the lower affinity T state when no or little Ca 2þ is bound, whereas it shifts towards the higher affinity R state when plenty of Ca 2þ is bound [27], which causes the cooperative effect. The identical 38 sets of 14 experimental traces as used above were fitted with the MWC model. Here also, both cooperative sets were considered to be identical and the cooperative part of the model was set up such that k on(R) , k on(T) , the apparent K d (K d(app) ) for the pairs, and the Hill number (n H ), could be fitted. This was achieved by adding a calculation step that determined K d(R) and K d(T) from the latter two parameters (see Protocol S1): As discussed above, L is dependent on the number of Ca 2þ ions bound to CR. To establish this dependence, we changed the forward and backward rate constants between the R and T states equally: We started with the same values for the (apparent) K d 's and n H as above with various combinations of k on starting values between 10 5 and 10 8 M À1 s À1 . As with the first model, occasionally a set of traces was fitted with values that deviated significantly from the general population, but again we found that only with one general set of constants, all 38 sets could be accurately fitted; the details are reported in Table 2. For comparison, we depicted the fit and the fit errors obtained with the MWC model, using the identical set of traces used for the new model (compare Figure 3A-3G, blue vs. red traces;  (Figures 4 and 5). For the independent site ( Figure 4A and 4B), results based on the MWC model (blue symbols) were essentially identical to the ones found with the new model (red symbols, compare also Tables 1 and 2 for the independent site V). Also, the properties of the apparent K d and the n H of the cooperative sites were similar between the two models ( Figure 5A, compare green and pink circles, and Figure 5E, compare blue and red symbols). This is a first indication that both models quantitatively describe the same process. Obviously, the detailed descriptions for the cooperative sites, applying either the new or the MWC model, deviate from one another, based on the differently modeled processes as described in Figure 2.

Comparison of the Kinetic Properties of CR Determined by Either the New or the MWC Model
Both models accurately fitted the data in a quantitative manner; based on the quality of the fits, they were indistinguishable, and technically can both be used to quantify the kinetic properties of Ca 2þ binding to CR. In particular, results for the independent site (V) are virtually identical (Figure 4). For the cooperative sites, both models describe binding sites that have a similar steady-state/affinity profile ( Figure 5A; apparent K d and Figure 5E; n H ). Our uncaging experiments were performed over a fairly narrow range of resting [Ca 2þ ] (2.0-5.3 lM), dictated by the constraints that most (.99%) of the DMn should be in the Ca 2þ -bound form to obtain valuable data. At lower resting [Ca 2þ ], the unbound DMn, present at higher concentrations than the free [CR], would rapidly rebind most of the released Ca 2þ [25], making CR's relative contribution to the [Ca 2þ ] decay small and difficult to distinguish. Because buffering kinetics depend on both the on-rates and the concentration of free buffer, the overall Ca 2þ binding speed to DMn would be much higher than that to CR, thus masking Ca 2þ binding to CR. With much less Ca 2þ -bound DMn, the changes in [Ca 2þ ] would be quite small and difficult to detect. Such technical constraints do not allow performing the experiments over the whole ''physiological'' range of [Ca 2þ ], e.g., from 10 nM to 100 lM. Thus, we could not exclude that the two models describe systems with different kinetic properties outside the boundaries of our experiments. This possibility was tested by examining the behavior of each model as a filtering system for Ca 2þ signals regulated by resting levels of [Ca 2þ ]. The filtering properties of both models were determined over a range of conditions covering the whole physiological range that CR is expected to encounter. CR (500 lM) was subjected to a wide frequency range (0.3 Hz to 10 kHz) of small (1 nM) sinusoidal perturbations of the [Ca 2þ ] at a wide range of starting [Ca 2þ ] (1 nM to 100 lM). The resulting Ca 2þ ''waves'' were close to sinusoidal. We used their amplitude as output of the filter to determine the transfer function of CR ( Figure 6). The attenuation of the sine wave is plotted as a function of frequency and resting [Ca 2þ ]. Both models ''filter'' Ca 2þ signals in a very similar way; the signals are less attenuated as the frequency gets higher (CR acts as a high-pass filter) or when CR becomes fully occupied as the starting [Ca 2þ ] gets higher. Remarkably, both models show a similar strong increase in attenuation at the lower frequencies, when the Ca 2þ concentration gets close to the apparent K d value of the cooperative sites, i.e., approximately 1.5 lM. The similarity of the transfer function of CR using either model indicates that they quantify the kinetics of Ca 2þ binding by CR in a similar way (but via different mechanisms) over the whole physiological range of conditions.

Comparing the New Model and the MWC Model
Both our new model, which is closely linked to the Adair-Klotz model [6,7], and the MWC model can be used equally well to quantify the Ca 2þ -binding kinetics of CR. However, we consider the new model to facilitate the ''intuitive'' understanding of how the kinetic properties of the cooperative sites relate to the binding kinetics at the level of the whole protein or at the macroscopic level. The Adair-Klotz model is the most general description of equilibria in terms of stochiometric binding. It describes the steady-state equilibrium using the constants (K 1 , K 2 . . ..K n ) for the successive binding (or macroscopic) steps, but not as the affinity constants of the individual (or microscopic) binding sites: for which in equilibrium, the fractional occupation (m) of a protein P is described by the Adair-Klotz equation: Usually, rate constants are not denoted in the macroscopic equilibrium equation (Equation 15); instead, only the K values are denoted, which is sufficient for steady-state descriptions. The equilibrium constants for the new model (K T and K R ) and the MWC model (K T , K R , and L 0 ) can often rather easily be translated into macroscopic K values [10] (also see Protocol S1). Therefore, it is fairly simple to relate any of the steadystate constants of cooperative models (new, MWC, and KNF) to the more generally used Adair-Klotz equation.
In addition to the calculation of the steady-state equilibrium (Equation 17), the macroscopic Adair-Klotz model compiles the binding of multiple binding sites into an intuitively easy-to-understand sequential binding model. It would even be more insightful if one could also obtain the rate constants for each macroscopic step. Unfortunately, the macroscopic rate constants are generally extremely hard to define when cooperative mechanisms are involved. For example, the macroscopic k on(1) for the MCW model depends on the relative amounts of totally unoccupied molecules in the R and T states (see Figure 2B). At steady state, this equilibrium is fairly straightforward, as this is simply defined by L 0 . However, when the balance is disturbed by a sudden change in Ca 2þ concentration, it will disturb the equilibrium between unoccupied molecules in the R and T states. This equilibrium will settle over some time according to k þ and k À . During this time, the relative amount of binding sites in states R and T is dynamically changing, making the macroscopic k on(1) itself dynamic. With most cooperative models, the macroscopic rate constants will be dynamic because they are dependent on most perturbations. This makes the rate constants very difficult to interpret (and to calculate).
However, for the new model of cooperativity, the macroscopic rate constants are easily defined and are truly constant. For instance, for two cooperative sites as described in this paper: and Through these simple relationships and according to our data, CR can be quantitatively described as a mixture of two ''virtual'' CaBPs. The cooperative part can be described as: and the independent part as: New Insights into the Regulation of Cellular Ca 2þ by CR As described above, at a starting [Ca 2þ ] around the apparent K d for the cooperative binding sites, CR will more effectively buffer perturbations at lower frequencies ( Figure  6). Thus, the Ca 2þ -buffering kinetics of CR clearly depends on the starting [Ca 2þ ] that determines the distribution between states T and R of the cooperative binding sites. While more cooperative sites get occupied, more Ca 2þ binding will take place through the second faster binding step as described in Equation 20. To better understand the cooperative nature of Ca 2þ binding by CR, we simulated with the new model a 1 lM step in [Ca 2þ ] from a resting [Ca 2þ ] of 10 nM in the presence of 100 lM CR. In comparison, we also simulated the widely used synthetic Ca 2þ buffers BAPTA (K d ¼ 160 nM, k on ¼ 2 3 10 8 M À1 s À1 , Maxchelator software version 10/02, see Materials and Methods) and EGTA (K d ¼ 70 nM, k on ¼ 1 3 10 7 M À1 s À1 [2]). Under these conditions, the [Ca 2þ ] decay kinetics mediated by 100 lM CR could be faithfully reproduced by either 7.8 lM BAPTA or 153 lM EGTA ( Figure 7A). Now, to compare the binding kinetics of CR to the two synthetic chelators without cooperative binding, we kept all parameters constant, i.e., 1 lM steps in [Ca 2þ ], buffer concentrations of CR, BAPTA, and EGTA, and only varied the resting [Ca 2þ ] (10 nM to 10 lM) from which the [Ca 2þ ] step was induced. As the resting [Ca 2þ ] increases, the [Ca 2þ ] decay kinetics in the presence of either one of the noncooperative chelators BAPTA and EGTA is slowed down (Figure 7A-7G). This slowing results from the fact that at higher initial [Ca 2þ ], more BAPTA or EGTA molecules will be in the Ca 2þ -bound form,  and thus, less Ca 2þ -free buffer is available. In contrast, as the initial [Ca 2þ ] increases from 10 nM to approximately 1.1 lM, CR speeds up the decay in [Ca 2þ ] and buffers Ca 2þ faster until approximately 1.1 lM resting [Ca 2þ ], above which CR behaves ''similarly'' to EGTA or BAPTA. More simulations have indicated that the exact breaking point between this novel behavior and classic behavior is dependent on the CR concentration and the Ca 2þ step size; however, the breaking point is always close to 1.1 lM for CR. Evidently neither EGTA nor BAPTA are able to mimic the properties of CR over the whole physiological range of [Ca 2þ ]. Yet for a given resting [Ca 2þ ] and a defined step increase in [Ca 2þ ], a BAPTA or EGTA concentration can be found that will closely mimic the effect of CR on the [Ca 2þ ] decay kinetics. But with only a slight change in the initial [Ca 2þ ] or the step size, this particular concentration of the synthetic chelator will not accurately reflect the action of CR. As an example, the concentrations of BAPTA or EGTA needed to mimic Ca 2þ binding by CR for the simulations (Figure 7A-7G) are shown in Figure 7H. Since CR has five binding sites, 100 lM CR is, in terms of Ca 2þ -binding sites, equivalent to 500 lM of either EGTA or BAPTA. For the example shown, at initial [Ca 2þ ] smaller than 0.3 lM, the concentration of EGTA needed to mimic 100 lM CR is in that same order of magnitude (;150-1,400 lM). The amount of BAPTA necessary to mimic CR is on the order of 500 lM when the initial [Ca 2þ ] is approximately 1-10 lM (;150-800 lM, respectively) ( Figure  7H). Thus, one can conclude that CR behaves EGTA-like around physiological resting [Ca 2þ ] (20-100 nM) typically seen in neurons and more BAPTA-like at the higher [Ca 2þ ] observed during bouts of activity.

Quantifying Kinetics of Cooperative Binding
The kinetic properties of Ca 2þ binding to CR were quantified using two different models featuring cooperativity between Ca 2þ -binding sites. The quantification of ligand binding to a protein with multiple binding sites is inherently difficult. First of all, when assuming that all binding sites are never truly identical, then for each binding site, at least two constants need to be determined. When these binding sites cannot be studied individually, the number of variables to be simultaneously determined quickly increases in a mathematical description of the protein. Consequently, degrees of freedom in the model will increase, which will decrease the accuracy and enhance the variability of the fitting. These difficulties are further exacerbated by cooperative binding, in which properties of the binding sites dynamically change, resulting in even more variables to determine. Although it is unlikely that any two binding sites are truly identical, a general way to decrease the number of variables is to assume the different binding sites to be identical in binding and allosteric behavior [12,13]. Even with such simplifications, the number of variables remains fairly large. To overcome this problem, we performed our experiments at a variety of initial conditions and at a variety of uncaging energies. In this way, we created a large set of data that varied by known parameters. We used a bootstrap method to create several datasets to be fitted simultaneously. The fit routine was set up in such a way that the models were required to consistently describe the properties of CR simultaneously over the whole range of variations in experimental conditions within one set, as well as over the whole range of uncaging energies. This technique resulted in (log)-normally distributed results so that erroneous fit results could be easily identified among the fits to any individual set.
Initial studies on CR have revealed the protein to contain cooperative Ca 2þ binding sites as evidenced by Hill coefficients of 1.3 or 2.4 [4] and 1.9 [22] for human and chick CR, respectively. Thus, this CaBP appeared to be well suited to serve as a model to quantitatively describe its cooperative Ca 2þ -binding properties. Even if such a description would fail to unravel the ''exact'' physical interpretation of Ca 2þ binding to CR, the obtained data were expected to serve as a powerful tool to understand the role of CR in shaping intracellular Ca 2þ dynamics in neurons. To quantify the cooperative binding to CR, two models were tested: a newly developed one and the MWC model. The latter was chosen because it has been used earlier to determine cooperative binding to biologically active molecules [14,16]. These studies were limited to special cases where transition rates between allosteric isoforms are much slower than the binding rates [14,15] or where binding and unbinding rates could be measured independently [16]. Here, we showed that it is feasible to determine the rate constants of CR using a MWC model outside of these constraints. In addition, we developed a different model closely related to the Adair-Klotz model [6,7], in which we did not account for the transition rate between the two possible states for a cooperative binding site. Both models yielded a similar quantitative description of CR's cooperative properties, but as long as crystal structures of CR in different states of Ca 2þ occupancy remain unknown, the exact physical interpretations of Ca 2þ binding to CR will not be available.

Properties of Cooperativity
Within one molecule, cooperativity can only be established by sets of at least two binding sites that change their properties based on the occupancy by Ca 2þ . In both models used here, four cooperative binding sites can occur in either the T state with low affinity for Ca 2þ or R state with high affinity for Ca 2þ . The dualistic nature of the binding sites causes the classical cooperativity as seen in steady-state binding studies in which the four cooperative sites in CR can be described with a K d of 2.5 lM and a Hill coefficient of approximately 2 [4,22]. The steady state properties for the cooperative sites of either the new model ( Figure 5A and Table 1, K d(app) ¼ 1.4 lM; and Figure 5E and Table 1, n H ¼ 1.9) or the MWC model ( Figure 5A and Table 2, K d(app) ¼ 1.5 lM; and Figure 5E and Table 2, n H ¼ 1.9) are in close agreement with the steady-state values for CR. Also for the independent site, the K d 's found with the new model ( Figure 4A and Table  1, 36 lM) or MWC model ( Figure 4A and Table 2, 41 lM) are in agreement with the earlier determined value (53 lM) by steady-state measurement [4]. The mechanism to create positive cooperativity with an initial low-affinity binding step as expressed in Equation 20 for the new model can be quite confusing. Since the first Ca 2þ -binding step is with a binding site in the low-affinity T state, it may appear that this will be limiting for the whole process, so that overall binding will only occur at higher [Ca 2þ ]. However, it should be considered that even at lower [Ca 2þ ] (e.g., at 100 nM inside a cell), at equilibrium. some of the sites in the T state will bind Ca 2þ , yielding some cooperative pairs with one Ca 2þ bound, which then rapidly leads to some fully occupied cooperative pairs: Since the second step in this process is governed by a highaffinity site, a considerable number of the cooperative pairs will be fully occupied. The combination of these reaction steps gives rise to the intermediate apparent K d (K d(app) ) (also see Equation 6 and, in Protocol S1, Equation S76). Although high-affinity R sites (short for ''a site in the R state'') only become available after Ca 2þ binding to CR's low-affinity T sites, they still play a determining role for the steady-state equilibrium. In dynamic situations, the slower binding of Ca 2þ to a T site has to precede binding of a second Ca 2þ to a faster R site, the former step apparently being rate limiting. However, at a given initial [Ca 2þ ], CR molecules are present in different states (metal-free, T state, R state, and completely Ca 2þ -bound) according to the parameters described in Tables 1 and 2; one example is described in more detail here. Assuming a 10 lM step in [Ca 2þ ] from a resting [Ca 2þ ] of 100 nM in the presence of 100 lM CR, the initial ratio for [unoccupied T sites]/[unoccupied R sites] is 395 lM/1.4 lM (calculated with Equations S34-S36 in Protocol S1). Of the 10 lM increase in Ca 2þ , almost all will be buffered by CR: 5.6 lM will bind to T sites, 3.9 lM will bind to R sites, and 0.35 lM to the independent site V. Evidently, the initially available R sites (1.4 lM) will not be sufficient, and most of the bindings to the R site will be time-limited by initial binding to T sites. However, at a resting [Ca 2þ ] of 1 lM, the initial ratio of unoccupied T/R sites is 250 lM/9 lM. In this case, 5 lM Ca 2þ will bind to T sites, 4.8 lM will bind to R sites, and 0.15 lM to site V. Thus, enough sites (9 lM) in the fast R state will be available and allow for fast buffering, not necessitating the slower stepwise change from the T to the R state. Model simulations of these experiments indicate that the surplus of sites in the R state will, within approximately 1 ms, bind up to 0.7 lM Ca 2þ more, which at a later time shifts to a CR molecule in the T state. By acting as a temporary substitute, the free sites in the R state can even ''speed up'' the eventual buffering of binding to sites in the T state. Albeit the initial amount of free sites in the R state is relatively small, it will significantly contribute to the buffering speed at initial [Ca 2þ ] around K d(app) . Also, the low-affinity independent site, which will hardly play a role in the steady-state buffering, initially binds up to 1.8 lM Ca 2þ which is later ''transferred'' to sites in the T state and will contribute to the overall buffering speed. However, this contribution is virtually identical when starting from either 100 nM or 1 lM Ca 2þ .

Unusual Nonlinear Properties of CR
To circumvent the difficulties of ''real'' CaBPs encountered in electrophysiological experiments (e.g., washout, unknown concentrations, and kinetic properties), they are often substituted by the artificial Ca 2þ buffers BAPTA and EGTA. Since BAPTA and EGTA have significantly different on-rates (2 3 10 8 M À1 s À1 and 1 3 10 7 M À1 s À1 , respectively), it is generally assumed that a process influenced by BAPTA, but not by EGTA, must have Ca 2þ -binding on-rates comparable to ones of BAPTA. So far, two studies [28,29] have inferred the kinetics of Ca 2þ binding by CR from comparing it to BAPTA and EGTA.
It was thus surmised that CR must have one or several binding sites with fast on-rate(s) [28]. The conclusion was drawn from the finding that addition of BAPTA, but not EGTA, to CR-deficient cells could rescue the CR deficiency. However, such generalizations are certainly error-prone, because both EGTA and BAPTA have K d values much lower than CR. This will considerably affect the results, since the speed of buffering is also dependent on the concentration of free sites and not only on the rate constants. Under resting conditions inside cells (;100 nM [Ca 2þ ]), at least 40% of either EGTA or BAPTA are occupied by Ca 2þ ions and do not add to the buffering speed. Our findings show that if BAPTA (or EGTA) can replicate a cellular buffering process under certain experimental conditions defined by, e.g., initial resting [Ca 2þ ], step size, and geometry of Ca 2þ influx, it does necessarily mean that the cellular buffering is ''comparable to BAPTA or EGTA'' and therefore ''fast'' or ''slow,'' respectively.
The exact intracellular distribution of CR is not well understood [30]; although principally considered as a cytosolic protein, a fraction of CR molecules could be anchored to specific sites [31][32][33], leading to higher local concentrations, possibly at Ca 2þ hotspots, as suggested for calmodulin around the L-type Ca 2þ channel [34]. Such local accumulation could lead to local high buffering speeds, especially at hotspots, where the higher local Ca 2þ concentrations would drive CR into a faster mode. Consequently, if CR is concentrated at certain subcellular compartments, the concentration of freely diffusing CR molecules will be lower. And if this freely moving CR is only confronted with smaller Ca 2þ signals, this will result in slower CR Ca 2þ -buffering kinetics. Slower properties of CR at lower [Ca 2þ ] together with its suggested mobility (a diffusion coefficient of approximately 20 lm 2 s À1 , assuming a similar mobility as the closely related CaBP CB [29,35]), supports the notion that part of CR will slowly bind and release Ca 2þ and thus act like a ''slow'' buffer. In contrast, the faster Ca 2þ buffering of CR at Ca 2þ hotspots or when present in the ''fast'' mode, i.e., with one of the cooperative sites in the Ca 2þ -bound form, may explain earlier findings that BAPTA could functionally rescue CR deficiency [28,36]. According to our findings, ''simple'' Ca 2þ chelators (EGTA and BAPTA) can never fully replicate certain functional aspects of CaBPs, because the complexity of Ca 2þ binding to ''real'' CaBPs such as CR cannot be mimicked by small synthetic Ca 2þ buffers lacking cooperativity.
The steady-state aspect of cooperative binding has been reported and analyzed in detail for Ca 2þ sensors such as calmodulin (for a review, see [3]). Cooperativity was also reported for CB [37,38], a protein with sensor and buffer functions [35,38], but the quantitative aspects of cooperativity have not yet been investigated in detailed steady-state binding studies. Nevertheless, cooperative binding has not been modeled in a CaBP to examine its effect on cellular Ca 2þ transients. Our results on CR pave the way to more realistically model intracellular Ca 2þ dynamics, thus leading to a better understanding of the spatial and temporal actions of Ca 2þ within a cell. The importance of correctly determining the physiological actions of CaBPs was recently shown in Xenopus oocytes. The effect of parvalbumin (PV), a CaBP with two Ca 2þ binding sites and minimal cooperativity [39], could be closely mimicked by the synthetic slow Ca 2þ -buffer EGTA [29]. Yet, the effect of CR on IP 3 -mediated Ca 2þ release was significantly different from that of the fast Ca 2þ -buffer BAPTA, in particular at low [IP 3 ], when Ca 2þ elevations were small. Under these conditions, CR caused a leftward shift in the concentration-response relationship as observed with the slow buffer PV. Under the same conditions, CR produced localized Ca 2þ transients or ''puffs,'' a phenomenon never observed in the presence of BAPTA [28]. Our findings strongly support the hypothesis [29] that the kinetic properties of individual CaBPs are finely tuned to specific cellular functions and may explain the need for a large number of CaBPs (more than 240 EF-hand-containing proteins) detected in the human genome [40]. Our novel approach determining the cooperative kinetics of Ca 2þ binding to CaBPs will lead to a better understanding of their highly specialized roles in cellular Ca 2þ signaling. In general, the method should also be applicable to any CaBP or to other multisite cooperative binding processes. This is expected to yield a more detailed understanding how CaPBs shape the spatiotemporal aspects of Ca 2þ signaling. Our findings may also help to reconcile some reported discrepancies concerning CR's function and putative effects on biological processes. The new model devised in the present study extends the analysis of cooperativity beyond the static steady-state condition, providing a powerful tool for the investigation of the dynamics and functional significance of cooperative binding in general. ratio of 40.0. These values were used in the mathematical model (see below) to describe the properties of the two batches of OGB-5N used in the various experiments. We have no explanation for the variability between these two batches other than the fact that specific contaminations might occur in different batches from the supplier (Molecular Probes, personal communication).
Properties of DMn. For each group of experiments, we determined the properties of DMn independently by uncaging experiments with no protein present, as described before [25]. These properties of DMn were then set for that specific experiment to compensate for possible differences between DMn batches. The observed properties of DMn were comparable to ones previously found (see Table 1 in Protocol S1) [25].
CR purification and determination of CR concentrations. Human recombinant CR was expressed in Escherichia coli and purified with a series of chromatographic steps as described before [4,44]. The purity of the isolated protein was estimated to be greater than 98% as judged from bands on SDS polyacrylamide gels (unpublished data). The initial protein concentration was determined by absorption measurements at 280 nm and using a molar extinction coefficient e 280nm of 26,860. Small aliquots of the protein (100-500 lg) in 10 mM (NH 4 )HCO 3 , 0.1 mM CaCl were lyophilized and then reconstituted directly in the solutions used for the uncaging experiments. To accurately determine the protein concentrations of all solutions used for the uncaging experiments, 10-15-ll samples were removed, stored at À20 8C, and simultaneously measured at the end of the series. The protein concentration was measured using a detergent-compatible assay based on a folin-phenol reagent (Bio-Rad) and using bovine serum albumin (BSA) as standard. All samples were measured in duplicates. Initial tests with solutions containing DMn and OGB-5N revealed that the colorimetric effect of these compounds was negligible at the concentrations present in the experimental solutions. The accuracy of the concentration measurements was validated by one round of fitting, in which the CR concentrations were fitted by the model, while the kinetic rates were allowed to deviate maximally 10% from their expected value. These fits confirmed the results of the protein assay.
Data analysis. All data were analyzed using MS Excel (Microsoft) and Berkeley Madonna 8.0 (University of California Berkeley). To determine the kinetic parameters (association and dissociation rates) from the fluorescence recordings, we used a mathematical model build in the ordinary differential equation solver Berkeley Madonna 8.0 that incorporates all of the reactions in the uncaging solution ( Figure 1). The DMn uncaging and OGB-5N signaling part of this model was used earlier to determine the exact properties of DMn [25]. This model was expanded with a part to simulate the binding of Ca 2þ to CR (see Figure 1B, and for a complete description of the models, see Protocol S1).
The model fit the simulations to the fluorescence recordings by iterating the parameters with the fourth-order Runge-Kutta method. The output of the model is: where F(t) is the fluorescence acquired over time (t) with t ¼ 0 at the time of the flash, and F rest is the resting fluorescence averaged over 50 ms before flash delivery . To more accurately fit the fast-rising phase while avoiding bias from late slow-decaying phase of the fluorescence transients, data points were omitted exponentially towards the end of every trace the fitting routine. Figure S1. Ca 2þ Transients Grouped According to the Experimental Conditions All individual Ca 2þ transients are grouped according to the experimental conditions as mentioned in the table. The red traces are the traces shown in Figure 3 of the paper. In the table for the OGB-5N column, I refers to lot number 34B1-2 and II to lot 15C1-2. Found at doi:10.1371/journal.pbio.0050311.sg001 (594 KB PDF). Figure S2. Pseudo-Random Picking of Fit Sets Flowchart of the compilation of 38 randomly selected sets of 14 traces derived from seven groups of data shown in Figure 3 of the paper. Random sets of traces were picked from the seven groups of data; from every group, two traces were picked per set. For the 38 sets, 76 picks are needed from every group. This means that every trace has to be picked 76/n times (n's are not equal for each group, see table), if every trace of a group is to be picked an equal number of times. Since 76/n is most likely not an integer, we picked every trace at least X times, where X is the closest smaller integer than 76/n (see table). To reach the number of 76 picks, 76 À n 3 X ¼ Y, see table) traces have to be picked one more time (X þ 1 times in total). By picking the traces in the way described here, we ensure that each trace within a group is used at approximately the same frequency for the fit sets.