Modulation of Calmodulin Lobes by Different Targets: An Allosteric Model with Hemiconcerted Conformational Transitions

Calmodulin is a calcium-binding protein ubiquitous in eukaryotic cells, involved in numerous calcium-regulated biological phenomena, such as synaptic plasticity, muscle contraction, cell cycle, and circadian rhythms. It exibits a characteristic dumbell shape, with two globular domains (N- and C-terminal lobe) joined by a linker region. Each lobe can take alternative conformations, affected by the binding of calcium and target proteins. Calmodulin displays considerable functional flexibility due to its capability to bind different targets, often in a tissue-specific fashion. In various specific physiological environments (e.g. skeletal muscle, neuron dendritic spines) several targets compete for the same calmodulin pool, regulating its availability and affinity for calcium. In this work, we sought to understand the general principles underlying calmodulin modulation by different target proteins, and to account for simultaneous effects of multiple competing targets, thus enabling a more realistic simulation of calmodulin-dependent pathways. We built a mechanistic allosteric model of calmodulin, based on an hemiconcerted framework: each calmodulin lobe can exist in two conformations in thermodynamic equilibrium, with different affinities for calcium and different affinities for each target. Each lobe was allowed to switch conformation on its own. The model was parameterised and validated against experimental data from the literature. In spite of its simplicity, a two-state allosteric model was able to satisfactorily represent several sets of experiments, in particular the binding of calcium on intact and truncated calmodulin and the effect of different skMLCK peptides on calmodulin’s saturation curve. The model can also be readily extended to include multiple targets. We show that some targets stabilise the low calcium affinity T state while others stabilise the high affinity R state. Most of the effects produced by calmodulin targets can be explained as modulation of a pre-existing dynamic equilibrium between different conformations of calmodulin’s lobes, in agreement with linkage theory and MWC-type models.


Introduction
Calmodulin is a ubiquitous calcium-binding protein involved in many cellular processes, from synaptic plasticity to muscular contraction, cell cycle regulation, and circadian rhythms. Structurally, it is organized in two highly homologous globular domains joined by a flexible linker [1]. Each domain contains two calcium-binding EF-hands that can undergo a transition between closed and open conformations, the latter favored by calcium binding [2,3]. The transition to the open state results in exposure of hydrophobic residues able to interact with numerous binding partners [4]. However, some targets interact preferably with the closed form [5].
The shift of calmodulin's calcium saturation curve in the presence of targets was studied experimentally by several groups [6][7][8][9][10]. Targets that markedly increase calmodulin's calcium affinity include the calcium-calmodulin dependent kinase II (CaMKII), protein phosphatase 2B (PP2B) and skeletal muscle myosin light chain kinase (skMLCK). It was shown that, in general, targets do not bind calmodulin exclusively in either calcium-saturated or calcium-free forms, but rather bind to both forms with different affinities. Numerous biologically relevant targets exhibit this behavior, such as skMLCK, CaMKII, and the NaV1.2 sodium channel [7,[10][11][12][13]. The binding domains that preferably interact with calcium-free calmodulin are called IQ motifs, a family of 14-residue sequences named after the two most frequent initial amino acids (usually isoleucine, followed by a highly conserved glutamate) [5].
Calmodulin's properties, such as cooperativity and affinity modulation by binding targets, can be explained by an MWC allosteric model. The first model of allosteric transitions was introduced in the seminal paper by Monod, Wyman and Changeux in 1965, which dealt with multimeric proteins whose subunits could undergo concerted conformational transitions, and also be modulated by target binding either to the R state, or the T state, in an exclusive fashion [14].The model was then further extended by Rubin and Changeux to describe molecules that could be modulated by binding partners capable of binding both conformational states, but with different affinites [15]. The hypothesis of the existence of distinct calmodulin conformations in thermodynamic equilibrium, which is crucial to an allosteric model, is supported by evidence of conformational dynamics of the protein in solution, with a time constant on the order of microseconds [16][17][18], and was also suggested by theoretical and computational analysis of coarse-grained Hamiltonians [19,20]. Structural studies also supported the hypothesis that a small fraction of calmodulin molecules can exist in a more compact conformation even in the presence of calcium [21].
Experimental studies with constrained mutants showed that the capability to switch conformation is necessary for cooperativity [3,22,23], which is compatible with an MWC model. Interestingly, conformational transitions are not common to all EF-hand based proteins, despite their high structural homology in calcium-free conditions. For example, calbindin-28k is structurally nearly identical to a calmodulin N-lobe, but its EF-hand remains closed after binding calcium [24]. It was proposed that the hydrophobic pockets of the open calmodulin lobes are stabilised by the unusual local abundance of the usually rare methionine residues (over 6% for calmodulin, against the 1% of the average protein) which have the highest flexibility, minimum steric hindrance and minimum solvation energy of all hydrophobic amino acids, and can therefore adapt to both segregated and and solvent-exposed conditions [4]. The flexible aliphatic side chains of methionine can also easily establish contacts with very diverse substrates, contributing to calmodulin's promiscuity [25]. Some examples of calmodulin's capability to bind very diverse targets is shown in Figure 1. Figure 1. Examples of calmodulin's capability to bind diverse targets. Calmodulin's N-lobe and C-lobe are in pink and light blue, respectively. Calcium ions are in green. Left: calcium-free calmodulin binding to a peptide of neurogranin (red) with its C-lobe [26]. Center: calcium saturated calmodulin wrapping around a peptide of MLCK (blue) with both lobes [27]. Right: calmodulin binding to a peptide of the SK channel (yellow) via a calcium-saturated N-lobe [28]. Several studies have shown that the linker region between the two domains is inherently flexible, and can assume different conformations depending on the orientation that the globular domains require for binding [29].
Despite their conceptual simplicity, the application of allosteric models can be challenging. The intrinsic calcium affinities of different conformational states are not easily measurable and need to be reverseengineered, because most experimental results are fitted using the phenomenological Hill or Adair-Klotz models, which do not incorporate conformational transitions. In addition, modelling intact calmodulin, within which each domain can switch its conformation independently, results in a large number chemical species that need to be explicitly enumerated (combinatorial explosion), and a larger number of parameters that need to be fitted simultaneously. A previous allosteric model of calmodulin was developed by Stefan et al. [30], to model the differential calmodulin-dependent activation of calcineurin and CaMKII in synaptic plasticity. However, the model postulated that both lobes would undergo concerted conformational transitions, and also had similar calcium-binding properties. In fact, the two lobes possess a remarkable degree of autonomy [31], and the calcium saturation curve of calmodulin was amost exactly recovered by superposition of the saturation curves of tryptic fragments TR1C and TR2C, containing respectively the N-terminal or C-terminal lobe only [6,32,33]. The four alpha-helices of each globular domains form two EF-hands that work together as one cooperative unit with two calcium-binding sites [34].
Despite of their high level of structural similarity, the two domains also exhibit strikingly different affini-ties and binding kinetics for calcium ions. The C-lobe has a 10-fold higher affinity, but much slower kinetics, than the N-lobe [35,36]. The clearly different saturation curves of the lobes, as observed in the intact molecule, are shown in Figure 2. . Different apparent calcium affinity of the N and C lobe. Experimental data was taken from references [3,9,22,35,[37][38][39] The differences in calcium-binding behavior are surprising in domains that are structurally so similar, but several studies showed that the EF-hand motif is tunable over a wide range of affinity and kinetics by modification of few key residues [40,41]. The incredibly high level of sequence conservation across species suggests that the different properties of the two lobes are some way crucial to calmodulin's function.
Moreover, even mutations that do not affect calcium-binding properties, but alter calmodulin's affinity for one or more targets, can be lethal [42]. The sequence of calmodulin's four EF-loops (12-residue calcium-binding pockets within each EF-hand) is given in Figure 3.
The C-lobe was reported to play a pivotal role in mediating calmodulin's calcium-dependent inter- Figure 3. Sequence of the four calcium-binding pockets (EF-loops) of Calmodulin. The calcium ion coordinates with the first,third,fifth,seventh and twelft residue of the loop, and also with the ninth through a coordinating water molecule (W). Highly conserved residues are in dark gray. . actions with its targets [43], and experiments with tryptic fragments of calmodulin consistently showed that the affinity of the C-lobe for calmodulin-binding domains is usually much higher than the affinity of the N-lobe [7,13]. Moreover, targets that bind calcium-free calmodulin seem to interact almost exclusively with the C-lobe [10,44,45]. In the light of these facts, we postulated that a major portion of the observed target-induced modulation effects could be explained by the interactions of the targets with the C-lobe only, and we focussed our initial effort on modelling the observed properties of TR2C, (i.e. the isolated C-lobe). Taking advantage of existing experimental data sets, we modelled the behavior of the TR2C tryptic fragment in the presence of peptides WFF and WF10 (peptides that express full-length and truncated versions of the calmodulin-binding domain of skMLCK) and Nav1.2IQp (a peptide based on the calmodulin-binding domain of the calcium-activated NaV1.2 sodium channel). We developed and parameterised a model of the tryptic fragment, in order to reliably reverse-engineer some intrinsic properties of the C-lobe of calmodulin. We then parameterised an MWC model of the isolated N-lobe by postulating that the R-states of both lobes had very similar affinities (as suggested by [6]), thus reducing the number of free parameters to estimate (see Methods).
This simplifying assumption was needed to circumvent the comparative scarcity of experimental data on the isolated N-lobe. The two submodels were then combined into a model of intact calmodulin, under the assumption that the protein behaved as the sum of its parts.
The MWC model explains cooperativity and target-induced affinity modulation as emergent prop-erties, rather than a priori assumptions, as opposed for example to the Adair-Klotz model. The two formulations are however mutually consistent, and MWC and Adair-Klotz models are interconvertible, as shown by Stefan et al. [46], since for any MWC model the corresponding Adair-Klotz parameters can be computed. A notable advantage of a MWC model is that the effects of multiple competing targets are straightforward to incorporate, simply by defining each target's affinity for the different conformational states of calmodulin. The work described in this paper shows that a carefully parameterised MWC model can indeed give reliable predictions of how individual targets modulate calmodulin affinity, and can also help investigate biologically relevant situations where numerous targets simultaneously modulate (and compete for) the same calmodulin pool.

Results
A diagram of the allosteric model of calmodulin with hemiconcerted conformational transitions is shown in Figure 4. Our model was built in three steps. First, we parameterised a model of TR2C, (isolated C-lobe), then parameterised a model of TR1C (isolated N-lobe), and finally the two submodels were merged into a model of intact calmodulin. (see Methods). The model of intact calmodulin (in SBML format) was deposited in BioModels Database [47] and assigned the identifier MODEL1405060000.

Preliminary analysis of possible parameterisations
The first requirement we set for our model was the capability to reproduce the saturation curve of TR2C alone. In the absence of targets, Equation 1 depends on only 3 free parameters (L, K R , K T ). We observed that the intrinsic affinity of the T state must lie within the wide, and biologically plausible, 1 µM -1 mM range. Once the value of K T is chosen, the model can be rewritten as a function of the "classical" MWC parameters, L and c, where c = K R /K T [14]. For any choice of K T , the score function of the fitting is a surface on the the (L, c) plane, with the same qualitative features, as shown in Figure 5.  Calmodulin has 4 calcium-binding sites (A, B, C, D) corresponding to one EF-hand each, and organised in two globular domains, N-lobe (pink) and C-lobe (light blue). Green arrows represent conformational transitions, red arrows represent calcium binding. Each binding site is modelled as a subunit that has access to two conformational states: open, high-affinity R state, and closed, low-affinity T state. Subunits on the same lobe are constrained to be in the same state, but the two lobes are allowed to be in different states. When no ligand is present, 4 possible configurations (TT, RT, TR, RR) coexist in thermodynamic equilibrium. A lobe that undergoes a T-to-R transition increases its affinity for calcium. Under these assumptions, it follows from classical linkage theory that ligand binding shifts the conformational equilibrium towards the open form, by stabilising the higher-affinity R state. Cooperativity is predicted as an emerging property.
of possible parameter choices that lie between two limit cases: in the first one, the change of affinity upon conformational transition is relatively small (c 0.01), and the allosteric constant is also comparatively small(L ≈ 100); in the second one, the change of affinity is much greater (c ≤ 0.001) and so is the allosteric constant (L ≥ 1000). In the first case, (higher c, lower L) the molecule has a high propensity to spontaneously switch conformation to the high-affinity R state, but the R state is stabilised less strongly Figure 5. Exploratory study of the general beahviour of a two-site MWC model. The affinity of the T state was fixed at specific values within the estimated range and the sum of the squared errors from the experimental points of the saturation curve, S, was plotted as function of the MWC parameters L and c. A lower score means a better fit. Over a very wide range of hypothetical, but physiologically plausible values of the T-state affinity (1mM − 10µM ), S(L,c) had the same qualitative behavior. All the pairs of (L,c) values that lie at the bottom of a flat valley (dark blue area) give a very good fit of the saturation curve of TR2C in the absence of targets. The knowledge of the target-free saturation function alone is therefore insufficient for univocal parameter identification.
by the binding of calcium. In the second case (smaller c, greater L), the opposite holds true: the molecule has a very low propensity to spontaneous conformational transitions, but the stabilising effect of calcium is much stronger. It must be stressed that in both cases the model is capable of giving an excellent fit of the experimental saturation curve in target-free conditions. The knowledge of the saturation curve in the absence of targets, alone, is therefore insufficient to discriminate between the two possibilities and univocally identify the parameters. The problem could be solved using the additional information provided by the saturation curves observed in the presence of targets, which allowed us to constrain the parameter space to be sampled during the fitting procedure. The obtained constraints implied that the first of the two limit scenarios mentioned above (with greater c and lesser L) had to be be discarded, because the resulting model would be unable to account for the extent of the target-induced affinity shifts observed in the experiments, as explained in the Methods section.

Truncation does not strongly affect the C-lobe's calcium affinity
We tested whether it was legitimate to reverse engineer the intrinsic properties of the C-lobe from those of TR2C. The level of saturation of the C-lobe can be monitored in both intact calmodulin and TR2C fragments by measuring the intrinsic fluorescence of their tyrosine residues [9,10,12,35,48], It was shown that truncation does not have a strong effect on the secondary structure and tyrosine fluorescence intensity of the C-lobe [49]. In order to verify that the calcium-binding properties of the intact calmodulin's C-lobe were maintained in the TR2C tryptic fragment, we compared the two saturation curves as presented in Figure 6. The data shows some scattering, and the uncertainty on the calcium concentration necessary to produce half-saturation is roughly of a factor two. This is probably due to slightly different experimental conditions (see the Discussion section). However all curves were similar and shared the same qualitative behavior, meaning that truncation does not dramatically alter the lobe's properties.
In the absence of targets, the COPASI model contains 3 independent parameters, (L,K R ,K T ). In the presence of target t, the saturation curve depends also on the target affinities K t r , K t T , which were both known quantities for WFF. The affinities of the peptide NaV1.2IQp for both the T and R sates were Figure 6. Comparison between the experimental saturation curves of TR2C and the C-lobe. Saturation of TR2C (red) and C-lobe of intact calmodulin (black), as a function of the free calcium concentration, monitored by intrinsic tyrosine fluorescence (data taken from references [7,9,12,22,32,35,49]). The available experimental data show some scattering, with an uncertainty of about a factor two for the concentration of free calcium necessary to produce half-saturation.
in the nanomolar range, and only estimated upper limits for their value were available in the literature (Table 1). We assumed that NaV1.2IQp affinity for the R state was equal to the estimated upper limit, while its affinity for the T state was set as an additional parameter to be fit, which we called K N aV In total, the 4 independent parameters K T , K R , K N aV T and L were fit simultaneously to the three experimental data sets.
The estimated parameters are summarised in Table 2, with the resulting saturation curves shown in Figure 7. The simulations of the model were overall in good agreement with experiments.
WF10 Truncated CaM-binding domain of skMLCK 1.1µM 712nM [7] NaV1.2IQp IQ domain of the NaV1.2 sodium channel ≤ 25nM ≤ 75nM [12] Table 1. Summary of the peptides used in the experiments described in references [7,9,12], and their affinities for the truncated C-lobe of CaM in the T and R state. The superscript t indicates calmodulin's affinity for a target. These peptides were used for the fitting and validation of the computational model.  [12] were measured with two different concentrations of peptide, and it was found that doubling the concentration of Nav1.2IQp (from 1:1.4 to 1:2.8 TR2C:peptide ratio) did not produce a further shift in the saturation curve, a behavior that our model was able to reproduce (see Supporting Information).A summary of all fitted and predicted saturation curves, plotted against the corresponding experimental data points, is given in Figure 8. Comparison between fitted and experimental saturation curves. Experimental data was taken from references [7] and [12].

The interplay of competing targets determines the affinity curve shift
The high calcium affinity of the C-lobe, further enhanced by interactions within some targets (CaMKII, PP2B, skMLCK), has led to speculation that calmodulin could activate such targets even at resting calcium concentrations, i.e. in the absence of calcium signals in neuronal or muscular cells [6]. We used the model of TR2C to perform a preliminary investigation of the effect of competing targets on calmodulin. At this stage of analysis we can only investigate the part of such interactions that are mediated by the C-lobe. However, as previously mentioned, the C-lobe was shown to be mainly reponsible for mediating calmodulin-target interactions, and targets that bind the T-state with high affinity, in particular those that appear to have very little interaction with the N-lobe [43]. We simulated the steady-state response of a system containing calmodulin and two competing allosteric targets using data taken from the literature, one binding prefentially to the the T state (T-state binding target, TBT) and Figure 8. Hill plot summarising the predicted saturation curves of TR2C in the presence of WFF, WF10 and NaV1.2IQp peptides, and comparison with the artificial datasets reverse-engineered from reference [7], and actual experimental data from reference [12]. In the range of physiologically plausible calcium concentrations (from nanomolar to tens of micromolar) the MWC model (continuous lines) is close to the Hill model (dash lines).
the other binding preferentially to the R state (R-state binding target, RBT). When both targets are present, the total saturation curve depends on their relative concentrations, as shown in Figure 9. Shifting the concentration of targets is a potential way to tune the saturation curve of the calmodulin pool.

Parameter estimation for the N-lobe
The model of the isolated N-lobe, in the absence of targets, was based on the model of C-lobe. The resulting parameters are summarised in Table 3. For the detailed procedure, see the Methods section. Each target alone can produce a marked shift in the saturation curve, but when both target are present in equal concentrations, the shift is much smaller. As targets, we chose peptides of two abundant neuronal proteins, neurogranin (Ng) and CaMKII, whose binding constants for TR2C were available in the literature. The chosen molar ratio (1:1:2.5) reflects estimated relative abundance of CaM, Ng and CaMKII in neuronal compartments. Their affinities for the T and R state of TR2C were respectively: 43nM and 1.05µM for neurogranin, and 88µM and 0.95µM for CaMKII ( [9,26]).

Model of intact calmodulin
The model of intact calmodulin was built by combining the two models of the N and C and lobe. The saturation curve agreed very well with the available experimental data as shown in Figure 10. A model including calmodulin and binding targets was then automatically generated in SBML format, as described in the Methods section, and simulated in COPASI.  Table 3. Summary of the estimated parameters for the isolated N-lobe. The R-state affinity was assumed to be equal to that of the C-lobe. The reported standard deviations were computed from the covariance matrix of the fit . Figure 10. Saturation curve of individual lobes and intact calmodulin, as predicted by our model, and comparison with experimental data. The datapoints for intact calmodulin were taken from references [6,7,50]. domain, corresponding to the portion that interacts with the C-lobe [7]. The saturation curve with WFF was markedly shifted to the left, while that in the presence of WF10 was biphasic.

Modulation of intact calmodulin by targets
Crucially, they were careful to measure the affinity of calmodulin and TR2C for these peptides, both in calcium and calcium-saturated conditions, which gave us an excellent initial estimate for the affinity of the peptides for the RR and TT state of intact calmodulin, and for the R and T states of TR2C. The affinity of each target for the isolated C-lobe in the R state was used as an estimate of the affinity of the target for the intact molecule in the TR state (i.e. we assumed that the target interacted with the C-lobe was much more strongly than with the N-lobe, in agreement with available experimental data).
The estimated affinities for the four states are summarised in Table 4. The agreement between simulation and experiment was very satisfactory ( Figure 11). Table 4. Estimated affinities of WFF and WF10 peptides for the different conformations of intact calmodulin, based on assumptions and actual affinity measurements from reference [7]. The affinity of WFF for the RR state was only given as upper limit (< 0.1nM ), and was assumed to be 0.05nM . The affinity of both peptides for the TR state was assumed to be equal to the peptide affinity for the isolated C-lobe in the R state. In the case of WF10, which is supposed to only interact with the C-lobe, the peptide affinity for the RR and TR states were set to be equal. Each peptide's affinity for the RT conformation (which is expected to be underpopulated in steady-state conditions and therefore play a marginal role in this analysis) was set equal to the affinity for the TT state. Figure 11. Effect of skMLCK peptides on the saturation curve of calmodulin, as predicted by our model, and comparision with experimental data by [7].

Calmodulin modulation by individual targets
We have shown that allosteric regulation can explain the modulation of the calcium-binding properties of the TR2C fragment (and hence, of the C-lobe of calmodulin) by several calmodulin-binding peptides.
The immediate consequence is that for a given concentration of target, the parameter that determines the extent of the saturation curve shift, at steady-state, is the ratio of the affinities of the target for the R and T states of calmodulin. Plotting the experimental points against the fitted curves in the Hill plane also offer an explanation as to why the MWC and Hill models can both provide a very good fit, despite predicting two quite different qualitative behaviors: in the Hill plane, the Hill curve is a straight line, whilst the MWC model predicts a gradual shift from a lower to a higher asymptote. As shown in Figure 8, in the range of physiologically plausible calcium concentrations (from nanomolar to tens of micromolar) the MWC model is close to the Hill model. The Hill coefficient predicted by the MWC model also decreases in the presence of targets, in qualitative agreement with experiments ( Table 5). The simulations of intact calmodulin in the presence of WF10 shows that biphasic saturation curves can be produced by targets that only stabilise the high-affinity form of one lobe. WF10 is an artificial peptide, but targets that induce differential modulation of calmodulin domains can be found in nature, as in the case of the anthrax edema factor [51].
An allosteric model of calmodulin was previously developed in our group that could satisfactorily reproduce the dose-response curve of wildtype calmodulin and also explain differerential activation of PP2B and CaMKII during synaptic tetanic stimulation [30].

Calmodulin modulation by competing targets
Calmodulin in vivo is always in the presence of a large number of targets, which simultaneously modulate its activity. We used our model to investigate the effect of mixtures of competing targets on calmodulin's C-lobe. We exploited the seemingly predominant role of the C-lobe in mediating calmodulin-target interactions to test the behavior of our model in the simultaneous presence of two different targets that had higher affinity for the T state and the R state, respectively. As an example scenario we chose peptides of abundant neuronal proteins, whose binding affinities for TR2C were available in the literature. We chose a 1:1:2.5 molar ratio for the three proteins, (a scenario where the concentration of calmodulinbinding proteins is much higher than the concentration of calmodulin [52]), to account for the higher concentration of targets that favour the R-state. In the chosen conditions, the effect of the competing targets was almost cancelled out and the resulting calcium saturation curve was close to that in the absence of targets, as shown in Figure 9. Regulating the relative abundance of targets can therefore bidirectionally tune the amount of calcium bound to calmodulin's C-lobe.

Two-state approximation and parameter estimation
A two-state model is of course a simplification of something as complex as a protein's conformational dynamics. It is more plausible that calmodulin is capable of sampling a wider ensemble of conformations, and its high conformational plasticity allows it to regulate downstream protein targets that are structurally very diverse [10,25,53]. Moreover, the PEP-19 peptide, expressed in the cerebellum, was shown to regulate mostly calmodulin's calcium binding kinetics, with little effect on affinity, and could do so even when calmodulin was already associated to CaMKII [54]. These observations imply that the tuning of calmodulin's affinity and kinetics is highly adaptable, but capturing every possible feature of this mechanism was outside the purpose of the present work. Moreover, for the conditions focussed on here, the agreement of a two-state model with experimental results was already satisfactory. On a technical level, parameter estimation is made challenging by the uncertainties and the scattering shown by the available experimental data. Ideally, allosteric models would have to be parametrised using data sets produced in highly standardised conditions, designed to be as close to the in vivo conditions as possible. However, different groups performed similar experiments and measured apparent calcium affinities that sometimes differ by up to a factor of two or more. This could be due to a number of reasons, including variable levels of purity of reagents or proteins, or slightly different buffer conditions. For example, increasing KCl concentration from 92 to 152mM was reported to decrease the apparent calcium affinity by a comparable factor [55]. As a consequence, pinpointing the exact value of some parameters can be challenging. For example, the standard deviation on the fitted value of the allosteric constant is rather large. On the other hand, this is mostly due to the low sensitivity of the saturation function with respect to this parameter: even an estimation error of a factor two (which is still relatively small, and close to noise level) would shift the saturation curve by an amount that is comparable to the scattering in the available literature data, as shown in Figure 12.

Conclusions
We have shown that a classic MWC allosteric model can successfully fit the saturation curves of calmodulin lobes and of the intact protein, and also account for the effects of several biologically relevant targets.
This wide applicability was also achieved with a comparative economy of hypotheses (independence of two domains, each in thermodynamic equilibrium between two conformations that have different affinities for the ligand and also different affinities for each target), thus providing a useful conceptual framework upon which further modelling work can be constructed. A crucial point is that the apparent affinity of calmodulin is modulated by simultaneous dynamic equilibria with different targets that exhibit a preference (tighter binding) to one of the possible conformations of calmodulin.

MWC model of an isolated calmodulin lobe modulated by targets
Each isolated lobe of calmodulin exhibits cooperativity between its two calcium binding sites. We chose to model the isolated C-lobe (TR2C tryptic fragment) as a two-subunit, two-state MWC molecule, where the subunits undergo concerted conformational transitions. As a starting point, we assumed that both calcium-binding sites on the C-lobe were identical, as was observed for example for the CaM85/112 mutant [22]. Therefore, in our model, both sites have calcium affinity K T when they are in the T state, and a higher affinity K R when the lobe is in the R state. Targets can bind both conformations of the molecule, but with different affinities, as shown by several experimental studies [7,12]. A diagram of the resulting model is given in Figure 13.
Under the simplifying assumption, used throughout this work, that the two calcium-binding sites are identical, the following equations can be derived for the saturation function of TR2C in the presence of an allosteric target A [15]: where: L = T 0 /R 0 is the allosteric constant in the absence of targets; L is the allosteric constant in the presence of targets; α = [X] f ree /K R is the ratio between the ligand affinity of the R state and the concentration of free ligand; c = K R /K T is the ration of the ligand affinities of the R and T state; is the ratio between the concentration of free allosteric target and the target affinity of the R state; e = K t R /K t T is the ratio between the target affinity of the R state and that of the T state.
The above equations allow the calculation of the saturation level of TR2C in the presence of a known concentrations of free ligand (calcium) and allosteric targets (peptides). They also clearly show that the effect of targets is equivalent to a modulation of the allosteric constant, whilst the other parameters remain unaffected. The functionȲ describes a surface in the plane (α, γ), as shown in Figure 14. The ligand saturation of an allosteric molecule, at equilibrium and at a given concentration of free ligand and free target, is entirely determined by the following parameters (note that all affinities are expressed as dissociation constants): Figure 13. Diagram of the model used to represent isolated lobes of calmodulin (TR1C and TR2C). Each lobe has two binding sites for calcium, two possible conformations, and can bind targets with different affinities depending on its conformation. As a simplifying assumption, the two calcium-binding sites were assumed to have identical affinities. The molecule can exist in two conformational states, T and R. In the absence of ligand, the equilibrium constant for the R to T transition is the allosteric constant L. The R state has affinity K R for calcium and K t R for target t. The T state has affinity K T for calcium and K t T for target t. In accordance to classic linkage theory, each calcium ion bound to the molecule shifts the conformational equilibrium towards the state with the higher calcium affinity, causing a scaling of the allosteric constant L by a factor c = K R /K T . In an analogous fashion, the binding of a target shifts the conformational equilibrium towards the state that binds t with the higher affinity, and scales L by a factor e = K t R /K t T . The model of intact calmodulin was obtained by combining two of these sub-models. .
-the allosteric constant L (or isomerisation constant), defined as the concentration ratio between the R and T states in the absence of ligand.
-the ligand affinity of the binding sites in the T state, K T ; -the affinity change upon transition from T to R state, c = K R /K T ; -the affinity for the target when the molecule is in the T state, K t T ; Figure 14. Saturation function of an MWC allosteric model of a two-state molecule with 2 binding sites, such as TR2C, in the presence of a target with preference for the R state. Increasing concentrations of target produce a left shift in the calcium saturation curve.
-the change of affinity for the target upon transition from the T state to the R state, e = K t R /K t T .
A diagram of an MWC model of TR2C is given in Figure 13. However, most available experimental data were obtained performing calcium titrations of calmodulin the presence of a fixed total concentration of target, and therefore Equation 1 was not directly applicable for fitting purposes. This apparent difficulty can be circumvented analytically, by calculating the corresponding concentrations of free target, as shown for example by Martinez et al. [56]. However, we found convenient to follow a computational route, directly implementing the model in the biochemical pathway simulator COPASI, which supports parameter fitting of steady-state properties when a system is initialised with total concentrations of reagents in a reaction compartment. The analytical model was instead used for a preliminary study on the general behavior of an MWC model for several choices of parameters, and also to discriminate between alternative possible parameterisations, as shown in the Results section.

Choice of experimental data for parameter fitting and validation
Several published steady-state titration curves described how calmodulin and TR2C fractional saturation at varying concentrations of free calcium change in the presence of different peptides that mimic the consensus motif of several in vivo binding partners of calmodulin. Published datasets were digitised using the freely available PlotDigitizer program. The datasets of Bayley et al. [7] and Shea et al. [12] were of particular interest because they reported the affinity for the peptides both in the absence and at saturating concentrations of calcium. In an allosteric perspective, these data provide estimates of the intrinsic affinity of the peptide for the T-state (predominant in calcium-free calmodulin) and Rstate (predominant in calcium-saturated calmodulin). These affinities were thus known quantities in our model. The properties of the peptides used in the above mentioned experiments are summarised in Table 1. It must be noted that none of the experimental datasets in the literature reported error bars or standard deviations for the measured datapoints, even when the datapoints were averages of measurements performed in triplicate [12]. In reality, all plotted datapoints are affected by uncertainty both on the x (free calcium concentration) and y axis (saturation) although free calcium concentration is usually known to high precision [6].
Fitting procedure and choice of constraints on the parameter space When choosing the initial conditions for the fitting procedure we assumed that our model must account for a wide range of observed behaviors. An initial exploratory study on different possible parameterizations was performed using the analytical MWC model described in the previous paragraphs, using custom Python scripts. These preliminary results were used to formulate additional constraints on the parametric space that the computational model was allowed to sample fitted to the experimental data describing the saturation curve of TR2C in the presence of targets. In the MWC framework, the apparent ligand affinity is produced by a mixture of two distinct populations of molecules, in the R-state and T-state.
Target peptides shift the saturation curve by differentially stabilising the R and T state. Therefore, any observed saturation curve, with or without targets, is bound to lie between two limit curves corresponding to populations of 100% R state, and 100% T state. The two fully-stabilised states also exhibit no cooperativity, because no further population shift is possible. These observations provide a simple way to determine lower bounds for K T and upper bounds for K R , which in turn provide an upper bound for c = K R /K T . The saturation curve for TR2C in the presence of Nav1.2IQp showed an apparent affinity of about 10µM and a Hill coefficient greater than 1, indicating that the T-state was not fully stabilised [12].
Therefore the affinity of the pure T-state must be lower than 10µM . On the other hand, data obtained with full-length calmodulin, showed that some peptides can increase affinity and decrease cooperativity, shifting the saturation curve to the left to an extent that implies that the calcium affinity of the pure R state must to be higher than 20nM [6], assuming that the C-lobe is responsible for the higher affinity when CaM is bound to a target, as shown by data published by Shea et al. [9]. Taken together, these observations led to the following constraints on the allowed parameter values: The chosen score function to minimise was the sum of squared residuals between the predictions of the model and the measured values. The fitting was performed using the built-in parameter estimation functions of COPASI, using 1000 iteration of a genetic algorithm with stochastic ranking [57], with a population size of 20, and with the constraints on the parameter values described in the previous paragraph. A genetic algorithm is a non-local optimisation method and is therefore less prone to converging to local minima in comparision to gradient-based methods.

Calculation of Hill coefficients
The Hill coefficients of calmodulin and calmodulin in the presence of targets were calculated as described in Ref. [58], as the slope of the saturation function in the Hill plane: where α = [Ca 2+ ] f ree /K R is the free ligand concentration normalised by the ligand affinity of the R state, and the fractional saturation function, Y , is defined as the ration between the number of occupied binding sites and the total number of binding sites, at a given concentration of free ligand, and can be calculated as: For an MWC model with two identical subunits, the Hill coefficient as a function of ligand concentration describes a symmetric bell-shaped curve, with its maximum at the point of half-saturation.

Parameterization of isolated N-lobe
The model of isolated N-lobe (corresponding to the TR1C tryptic fragment) was formally identical to the TR2C fragment described in the previous paragraphs. However, the parameter estimation required a different approach. We could not directly fit the model to saturation curves in the absence and presence of targets, because the affinity of the isolated N-lobe for targets is very low [7], and the resulting shift of the saturation curve very weak [12]. Moreover, targets that bind calmodulin in calcium-free conditions have little or no interaction with the N-lobe, and their effect on its saturation curve was negligible. First we assumed, for the sake of simplicity, that the two calcium-binding sites of the N-lobe are identical (as we did for the C-lobe). Based on the experimental evidence [6], that the two lobes have very similar affinity when bound to high-affinity targets, we made the additional assumption that, when in the R state, both lobes had the same calcium affinity. The affinity of the C-lobe in the R state was already known from the parameterization of the TR2C model. We estimated the calcium affinity of the N-lobe in the T state by fitting the saturation curve by Grabarek and coworkers for the NCaM41/75 mutant (where the N-lobe was constrained in a closed conformation by a disulfide bond) with a non-cooperative model with two identical binding sites. Knowing the affinity of both the R and T state, the allosteric parameter c = K R /K T was readily calculated. The only remaining free parameter in the model was the allosteric constant L, which was estimated by fitting and analytical MWC model onto the available experimental points for the saturation curve of the N-lobe [3,9,35].

Generation of the model in SBML format
The model of intact calmodulin was built by assuming that each lobe still behaved as it would in its isolated form. When referring to the conformation of intact calmodulin, we adopt the following convention: "RT" means that the N lobe is in the R state, and the C-lobe is in the T state. Since the two lobes have different affinities and therefore saturate at different calcium concentrations, asymmetric states are biologically plausible. However, the resulting kinetic model is much more complex than the submodels comprising only one lobe. Two largely independent lobes imply that the conformations are four: RR,RT,TR,TT. The binding sites are also 4 (A,B on the N lobe and C,D on the C-lobe), which leads to 64 possible states for calmodulin only (see Supplementary Material). With the addition of targets, combinatorial explosion implies that the number of equations in the kinetic model quickly increases to more than one thousand. Manual modification of such models and inclusion of additional targets would be error-prone and labour-intensive. On the other hand, the interaction rules underlying the model are simple, and the complexity is purely combinatorial, which means the model lends itself to be generated iteratively. Taking advantage of the existing libSBML library, we wrote custom Python scripts to automatically generate allosteric models of calmodulin in the presence of targets, in the SBML format. SBML is natively supported by COPASI, thus making the setup process of different kinetic models both faster and more reliable.

Affinity of targets for asymmetric conformations of calmodulin
Our model allows for hemiconcerted conformational transitions (i.e. the two lobes can have different conformations, but the EF-hands on the same lobe are always in the same state). This poses no problem for calcium binding events, because calcium only binds to individual sites, whose affinity is determined by their state. Targets, on the other hand, being bigger molecules (either peptides of whole proteins) bind to the whole molecule and usually interact with both lobes at once. Their affinity for a given conformational state of calmodulin can be expected to be an interplay of their affinity for the two lobes. A modelling challenge is posed by the fact that in the available experimental data, the affinity of calmodulin for a given target was usually measured only in two limit cases, in the absence of calcium (when almost all calmodulin will be in the TT state) and at saturating concentrations of calcium (when calmodulin will be mostly in the RR state). In our model we also needed the affinities for the two asymmetric states RT and TR. To overcome this limitations, we observed that targets that bind preferentially the calcium-saturated calmodulin exhibit a much stronger interaction with the C-lobe than with the N-lobe [9,49] and that the N-lobe alone can only bind very weakly to targets, and the affinity is in the millimolar range in the absence of calcium [7,12]. Collectively, these observations led to the following simplifying assumptions: -The affinity of the target for the RR state is equal to the affinity for CaM, measured at saturating calcium concentrations; -The affinity of the target for the TT state is equal to the affinity for CaM, measured in the absence of calcium.
-The affinity of the target for the TR state is equal to the affinity for the TR2C fragment (isolated C lobe) at saturating calcium concentrations.
-The affinity of the target for the RT state is roughly equal to the affinity for the TT state With this set of assumptions, the model could reproduce several experimental saturation curves (see

Results).
Supplementary material TR2C model validation Figure 15. Prediction of the behavior of TR2C in the presence of the NaV1.2IQp peptide, and comparison with experimental data not included in the original fitting procedure [12].
Parameters for the TR2C model where we only examined steady-state properties. The model was thus written using arbitrary kinetic constants, chosen to ensure that the ratio between the forward and backward rate constants of each reversible reaction was equal to the equilibrium constant calculated by recursive application of mass action law. The model for TR1C is formally identical to that of TR2C and is omitted. Figure 16. Prediction of the behavior of TR2C in the presence of the WFF peptide, and comparison with experimental data not included in the original fitting procedure [7].

Parameters for the model of intact calmodulin
The generic names "tbp" and "rbp" were used to desctibe peptides that exhibit a tighter binding for the T and R state, respectively.  [7,12], this paper Table 6. Summary of the model parameters. The kinetic rates for the conformational transitions were calculated based on the value of the allosteric constant estimated in this work, and the value of the exchange constant for the conformational fluctuations k ex ∼ 20000s −1 , estimated by Malmendal and coworkers for the C-lobe [16]. The exact value does not affect steady state properties, but a solid estimate will be importat for future work where kinetic properties will also taken into account.