Modelling intracellular competition for calcium: kinetic and thermodynamic control of different molecular modes of signal decoding

Frequently, a common chemical entity triggers opposite cellular processes, which implies that the components of signalling networks must detect signals not only through their chemical natures, but also through their dynamic properties. To gain insights on the mechanisms of discrimination of the dynamic properties of cellular signals, we developed a computational stochastic model and investigated how three calcium ion (Ca2+)-dependent enzymes (adenylyl cyclase (AC), phosphodiesterase 1 (PDE1), and calcineurin (CaN)) differentially detect Ca2+ transients in a hippocampal dendritic spine. The balance among AC, PDE1 and CaN might determine the occurrence of opposite Ca2+-induced forms of synaptic plasticity, long-term potentiation (LTP) and long-term depression (LTD). CaN is essential for LTD. AC and PDE1 regulate, indirectly, protein kinase A, which counteracts CaN during LTP. Stimulations of AC, PDE1 and CaN with artificial and physiological Ca2+ signals demonstrated that AC and CaN have Ca2+ requirements modulated dynamically by different properties of the signals used to stimulate them, because their interactions with Ca2+ often occur under kinetic control. Contrarily, PDE1 responds to the immediate amplitude of different Ca2+ transients and usually with the same Ca2+ requirements observed under steady state. Therefore, AC, PDE1 and CaN decode different dynamic properties of Ca2+ signals.


NMDARs
NMDARs are glutamatergic ionotropic receptors that mediate the influx of Ca 2+ to the cytosol of the hippocampal spine during the induction of long-term forms of synaptic plasticity 12,13 . Structurally, hippocampal neuronal NMDARs are composed by two GluN2 subunits, which are the subunits that bind glutamate, combined with two GluN1 subunits that bind to the co-agonists L-glycine and D-serine 14 . For simplicity, we assumed that NMDAR co-agonists are constantly present in the synaptic cleft, which is justified by the continuous release and uptake of L-glycine and D-serine from the extracellular environment 14 .
Therefore, they were not simulated explicitly and the gating of NMDARs in the computational model was controlled solely by the glutamate concentration. According to the native composition of synaptic NMDARs 12 , two types of receptors with distinct kinetic properties were included in the model, di-heteromerics receptors GluN1/GluN2A, which account for approximately 65% of the NMDARs in the hippocampal CA1 of adult animals, and the receptors GluN1/GluN2B that was considered to represent the remaining 35% of the total NMDARs 15,16 . A total of 20 NMDARs were included in the model 17 . Both receptor types were simulated with the same set of reactions, but using different rate constants for the state transitions 11 .
A unique feature of NMDARs is their mechanism of activation that requires both the binding of glutamate and the depolarization of the postsynaptic cell membrane to promote the release of magnesium ions (Mg 2+ ) that block the pore of their channels in a membrane potential (Vm)-dependent manner 12,18 . The depolarization-dependent removal of Mg 2+ occurs only from the open channel 18 . In consequence, NMDARs detect the coincidence between the presynaptic activation, which releases glutamate, and the postsynaptic membrane depolarization 13 . The kinetic model used to simulate both types of NMDARs in this work was based on a previous model developed to simulate the gating of NMDAR channels in absence of Mg 2+ 19 , which was recently expanded to incorporate the interaction of NMDARs with Mg 2+ (Supplementary Figure S11A) 9,10 . All the rate constants used to simulate each type of NMDAR were taken from published models [9][10][11] . The complete description of the reactions and parameters used to simulate the NMDARs is listed in Supplementary Table I (Reac1-Reac17).

Ca 2+ Dynamics
The mechanisms of Ca 2+ dynamics implemented in our stochastic model consisted of Ca 2+ influx, buffering and extrusion from the cytosol. We have not included Ca 2+ diffusion across the spine neck to its parental dendrite in our model because it accounts for less than 10% of the Ca 2+ clearance 20 .
The influx of Ca 2+ in the model occurs exclusively through NMDAR channels, and was given by the follow equation: (1) where NMDARopen stands for the number of open NMDAR channels unblocked by Mg 2+ (Supplementary Figure S11A), INMDA is the maximum unitary channel current (5 pA) 21 where Vm is the membrane potential, Erev is the NMDAR reversal potential (0 mV). However, since BioNetGen 26 , the software used to implement the model, does not allow the explicit simulation of equations such as (1), we simulated the influx of Ca 2+ in the BioNetGen as a first order reaction: The rate constant kinflux was estimated with a multiplicative factor (8.10 6 mol -1 .L) applied to equation (1)

PDE1
Eleven subfamilies of PDEs have been identified and, in all organisms, multiple isoforms of PDEs appear to be involved in the tight control of cAMP concentration 56 .
PDE1A2, a member of the PDE1 subfamily, is highly expressed in the brain 57 . PDE1A2 (referred here simply as PDE1) can interact with CaM fully or partially saturated with Ca 2+ , but it requires CaM bound to three or four ions to be stimulated 58,59 . CaM fully saturated with Ca 2+ interacts with PDE1 with affinity (KD) of 1 nmol.L -1 57 , and a Ca 2+ requirement of approximately 300 nmol.L -1 60 . The interaction between PDE1 and Ca 2+ /CaM was simulated with simple reactions of association/dissociation (Supplementary Figure S11B). The rate constants for these reactions were estimated from the KD and based on a slow half-life for the dissociation reaction, which was used to estimate a rate constant of dissociation of (kb) saturated with Ca 2+ was set as 1 µmol -1 .L.s -1 . This rate was kept unchanged for the interaction between PDE1 and CaM partially loaded with Ca 2+ . To simulate the interaction between PDE1 and CaM partially loaded with Ca 2+ , we assumed that the absence of Ca 2+ bound to the   (the terms I, II, III and IV   indicate the Ca 2+ -binding sites I, II, III Table   I: Reac86-Reac92).

cAMP dynamics and PKA
The balance between PDE1 and AC activity is fundamental to control the intracellular level of cAMP, one of the most important intracellular second messengers. The main target of cAMP is the cAMP-dependent protein kinase (PKA) 63,64 , an enzyme that frequently opposes CaN action. AC produces cAMP with a relative high basal rate 65  . This holoenzyme complex kept PKA C subunits in an inactive state in absence of cAMP, which binds to cyclic nucleotide-binding domains (CNB) located in PKA R subunits 64,70 .
The binding of cAMP to the two CNBs in each R subunit promotes a conformational change that leads to the release of the two active C subunits from PKA holoenzyme 64,70 . The R subunits occur in two distinct classes (RI, RII), each having an α and β isoforms 6,7 . RIIβ is the predominant isoform in the brain and the tetramer (RIIβ)2C2 has the highest cAMP requirement in comparison to other PKA holoenzymes 6,7 .
The crystalized structure of the full-length (RIIβ)2C2 solved recently indicates that each RIIβ subunit of the tetrameric PKA contacts the neighbouring R:C complex 6 . The physical contacts among different subunits in the tetrameric holoenzyme are responsible for PKA allosteric activation 6 . Several groups have developed computational models to simulate PKA activity 71,72 , but most of these models simulated PKA as a dimer of heterodimers, which is not supported by the experimental findings that clearly demonstrate that there is communication between the subunits of the tetrameric PKA during its activation 6,73 . To simulate PKA, we assumed that it is structurally composed by a dimer of two R subunits (RIIβ)2 coupled to two C subunits. The binding/unbinding of each C subunit to a RIIβ subunit of the dimer (RIIβ)2 was implemented as independent events in the absence of cAMP.
Endogenous PKA is autophosphorylated in the cells, which affects the affinities between the C and RIIβ subunits 10 . The phosphorylation of the RIIβ subunits by the C subunits was modelled as a first order reaction. PKA is largely in its autophosphorylated state in our model. To simulate the binding/unbinding of each C subunit to the phosphorylated RIIβ subunits, we used the same reactions implemented for the non-phosphorylated PKA, but different rate constants 10 .
cAMP can access both the CNB-A and CNB-B of the RIIβ subunits in the apoenzyme 7 . In consequence, we simulated the interaction of cAMP to each cAMP-binding site as random and independent events. The affinities and the rate constants observed for the binding/unbinding of cAMP to the isolated (RIIβ)2 dimer are different from the parameters reported for the full holoenzyme 74-76 . We assumed that each CNB of the tetrameric PKA has a 3-fold lower affinity for cAMP than the CNBs from the isolated (RIIβ)2 dimer 7,76 . The phosphorylation of the RIIβ subunits had no effects on their affinities for cAMP.
The binding of two cAMP molecules to each R subunit of the inactive PKA induces a decrease in the affinity between the C and RIIβ subunits, which leads to the dissociation of the holoenzyme into two free C subunits and the (RIIβ)2 dimer 64,70 . The KD for the interaction of the C and RIIβ subunit is 0.14 nmol.L -1 in absence of cAMP, and approximately 15-fold greater in presence of four molecules of cAMP 7 , which we implemented by changing the rate constant for the complex (RIIβ)2C2 dissociation. The activation of PKA fully saturated with cAMP is highly cooperative (nHill = 1.7). To capture this property, we assumed that, after the release of the first C subunit from the holoenzyme (RIIβ)2C2, the second C subunit was released with a faster rate constant. This cooperativity was implemented only for PKA fully saturated with cAMP.
Mutations of the CNB-B domain promote a decrease in the K1/2 for PKA activation, which is due mainly to a change in the rate constant (kb) for the dissociation of the complex (RIIβ)2C2 6,7 . We used this information to define the rate constants for the release of the C subunits from the tetrameric PKA in the absence of cAMP bound to the CNB-B sites: we kept kf for the complex association unchanged and multiplied the kb by twenty 6,7 . In contrast, mutations in the CNB-A domains promotes a strong increase in the K1/2 for PKA activation 7 , which we took into account to implement the dissociation/association between the C and the R subunits in absence of cAMP bound to the CNB-A sites. The rate constant for this dissociation was defined as the kb for the release of the C subunits from the tetramer (RIIβ)2C2 divided by fifth 7 , kf was kept unchanged. We simulated the interactions between the C and the phosphorylated RIIβ subunits using the same chemical reactions implemented for the non-phosphorylated RIIβ subunits, but different rate constants 10 .
To validate the thermodynamics properties of the PKA model implemented, we fitted its dose-response curves (equation (1)   The simulations began with the the total amount of subs completely phosphorylated.
Then, we stimulated the model using 100 glutamate pulses released at 1 Hz and 10 Hz to verified how the combined activity of CaN, AC, PDE1 and PKA regulates the dephosphorylation of subs (Fig. 5G-H).

Data Analysis
We