Analog Signaling With the “Digital” Molecular Switch CaMKII

Molecular switches, such as the protein kinase CaMKII, play a fundamental role in cell signaling by decoding inputs into either high or low states of activity; because the high activation state can be turned on and persist after the input ceases, these switches have earned a reputation as “digital.” Although this on/off, binary perspective has been valuable for understanding long timescale synaptic plasticity, accumulating experimental evidence suggests that the CaMKII switch can also control plasticity on short timescales. To investigate this idea further, a non-autonomous, nonlinear ordinary differential equation, representative of a general bistable molecular switch, is analyzed. The results suggest that switch activity in regions surrounding either the high- or low-stable states of activation could act as a reliable analog signal, whose short timescale fluctuations relative to equilibrium track instantaneous input frequency. The model makes intriguing predictions and is validated against previous work demonstrating its suitability as a minimal representation of switch dynamics; in combination with existing experimental evidence, the theory suggests a multiplexed encoding of instantaneous frequency information over short timescales, with integration of total activity over longer timescales.

Introduction

Many cellular inputs lead to transient changes in cytosolic calcium (Ca²⁺) levels, generating temporally complex signals that reflect a wealth of information (Berridge et al., 2003). As such, cells express highly conserved molecular decoders capable of translating Ca²⁺ oscillations into downstream signaling events that affect diverse processes such as gene transcription, development and aging, neural network homeostasis and the synaptic plasticity that underlies learning and memory (Lisman et al., 2002; Thomas and Huganir, 2004; Wen et al., 2004; Clapham, 2007; O'Leary et al., 2013; Tao et al., 2013; de Jong and Fioravante, 2014; Smedler and Uhlen, 2014). A celebrated example of a Ca²⁺ decoder is the protein kinase Ca²⁺/calmodulin (CaM)-dependent protein kinase II (CaMKII; Box 1), which can be driven by transient levels of cytosolic Ca²⁺ into either high or low states of switch-like activity. When stabilized through negative regulation by protein phosphatases, self-exciting (autophosphorylating) kinases such as CaMKII are an ideal component of signal amplification and have been previously likened to transistors on a computer chip, in that they may be turned on or off, presenting an ideal substrate for computation in cellular systems (Hunter, 1987; Ferrell and Ha, 2014).

Box 1. The bistable molecular switch camkii and synaptic plasticity.

Accounting for approximately 1–2% of all brain protein, CaMKII is a central hub of cell signaling networks and can exert both pre- and post-synaptic control over information transmission in the central nervous system (Lisman et al., 2002). Once bound to the Ca²⁺-CaM complex, the kinase's ability to cooperatively autophosphorylate can produce two distinct stable states: either high or low levels of enzymatic activation. Postsynaptically, after repetitive stimulation, the high activation state may persist after the Ca²⁺ signal subsides and can strengthen the connection between neurons, for example, the hippocampal CA3-CA1 synapses that support learning and memory (Lisman et al., 2012). However, it should be noted that the role of CaMKII autophosphorylation and constitutive activation is not fully understood or accepted (Michalski, 2013) and we are just beginning to gain better insight into the problem (Chang et al., 2017). This paper proposes that CaMKII's principal role is to meaningfully transmit information via its short term dynamics rather than store it permanently within levels of autonomously activated switch.

Presynaptically, CaMKII also modifies connection strength (Wang and Maler, 1998; Ninan and Arancio, 2004; Pang et al., 2010). In weakly electric fish, the αCaMKII isoform produces presynaptic potentiation in a motion sensitive, excitatory sensory feedback pathway (Wang and Maler, 1998; Clarke and Maler, 2017). The kinase also potentiates hippocampal CA3-CA1 synapses, as evidenced by knocking-out αCaMKII, which leads to reduced synaptic potentiation under paired pulse facilitation protocols when compared to the wild-type (Chapman et al., 1995). Through enzymatic phosphorylation of voltage gated Ca²⁺ channels and ryanodine receptors, αCaMKII can enhance Ca²⁺ entry and Ca²⁺-induced Ca²⁺ release in response to high frequency signals, potentially supporting hysteresis (Figure 1) and driving synaptic release (Catterall and Few, 2008). However, at the same CA3-CA1 synapses, post-tetanic potentiation protocols generate enhanced levels of potentiation in the same knock-out mice, illustrating that αCaMKII may also depress synaptic strength depending on the frequency and duration of the input (Chapman et al., 1995). Furthermore, αCaMKII has been shown to serve as a negative, activity-dependent regulator of neurotransmitter release probability at CA3-CA1 synapses (Hinds et al., 2003). This effect may be partially explained by the fact that CaMKII phosphorylates Ca²⁺-activated potassium channels that hyperpolarize the presynaptic terminal (Wang, 2008), decreasing the likelihood of Ca²⁺ entry and evoked neurotransmitter release. Intriguingly, αCaMKII also plays a non-enzymatic role in presynaptic CA3-CA1 plasticity by regulating the number of docked synaptic vesicles containing neurotransmitter (Hojjati et al., 2007). In this case, decreased transmitter release could be explained by the fact that αCaMKII is acting as a sink for intracellular Ca²⁺, lowering the cytosolic levels that drive the machinery of synaptic vesicle fusion and influencing the size of the readily releasable vesicle pool (Thanawala and Regehr, 2013; Jackman et al., 2016). The size of the readily releasable pool is directly correlated with release probability at hippocampal synapses (Dobrunz and Stevens, 1997), supporting a putative role for αCaMKII in control of presynaptic plasticity parameters via Ca²⁺ and CaM buffering (Hinds et al., 2003).

One of the most influential discoveries about CaMKII is its ability to decode the frequency of periodic Ca²⁺ pulses into distinct amounts of long lasting, autonomously activated kinase (De Koninck and Schulman, 1998). However, the interpretation of CaMKII as a frequency decoder has been criticized based on the fact that mean values of activity, evoked by different combinations of Ca²⁺ pulse size, duration and frequency, are ambiguously mapped into the same level of autonomously activated switch (Pinto et al., 2012), which suggests that the switch is actually integrating the Ca²⁺ input over longer timescales. Alternatively, this article focuses on whether the concentration of activated switch acts as a reliable (analog) signal that reliably encodes frequency information over short timescales (sub-seconds), where Ca²⁺ pulse size and duration are far more stable (Tank et al., 1995). The experimental evidence discussed above suggests that frequency coding by these “digital” molecular switches is more sophisticated than previously thought and that fast fluctuations in presynaptic αCaMKII around either the stable high- or low-activation state can better represent instantaneous frequency information, and, hypothetically, translate it into bidirectional control of synaptic strength in real-time.

FIGURE 1

Figure 1. Activation states of the bistable molecular switch model. (A) The model's potential function, U(x), visually describes the tendency for solutions to settle around one of two equilibrium points (x*), where the rate of change of switch activation, f(x), is 0 (parameters, r = 0.52 and c = 0.04). To the left of the stable equilibria (black circles), f(x)>0 (green), and to the right, f(x) < 0 (blue), which forces perturbations to settle back into those states. Conversely, the sign of f(x) is reversed on both sides of the unstable equilibrium (red circle), such that tiny perturbations push the switch away, toward either stable state. (B) As r or c change, f(x) changes and can result in the loss of bistability. (i) To illustrate, r is fixed as the input c is varied: small values only support low activation, but, as c grows, bistability emerges and eventually disappears as only the high activation state is supported when c>c_c (rightmost). A defining feature of bistability is the hysteresis effect, where the same value of a parameter may evoke different states depending on the history of activity. For example, the high activation state still exists for c less than the rightmost c_c and can only be lost when c falls below the leftmost c_c value. (ii) c is fixed while the negative regulation parameter r is varied. For small r, only the high activation state exists. As r grows larger, the system becomes bistable and, eventually, only the low state exists after crossing r_c. Panel (iii) shows a parametric plot of the critical values c_c(x) and r_c(x), that partition the parameter space, and the bifurcation surface summarize the analysis completely (iv).

The classic CaMKII experiments of De Koninck and Schulman provided the first demonstration that a molecular switch can decode the frequency of periodic Ca²⁺ pulses into distinct, persistent levels of high enzymatic activation (De Koninck and Schulman, 1998). Although experimental evidence still largely lacks for whether this persistent activation occurs within functioning cells (Michalski, 2013), there are recent indications that it does occur to some extent (Michalski, 2014; Urakubo et al., 2014; Rossetti et al., 2017) and that autophosphorylation is key to this process (Chang et al., 2017; Rossetti et al., 2017). Many modeling studies of CaMKII autophosphorylation dynamics capture the ability of the high activation state to persist beyond the original Ca²⁺ signal (known as hysteresis), which could potentially act over long timescales (seconds, minutes, and longer) (Zhabotinsky, 2000; Dupont et al., 2003; Graupner and Brunel, 2007). In these studies, the relationship between Ca²⁺ concentration and the state of the molecular switch are determined from simulations of detailed, parameterized systems of differential equations that are not readily amenable to deeper mathematical analysis; furthermore, these studies are restricted to periodic inputs and concerned with long timescale activation. In order tobetter understand frequency coding over short timescales (milliseconds to seconds) and its putative effect on synaptic plasticity (Box 1), this article analyzes a reduced description of molecular switch behavior when subject to noisy, aperiodic forcing, while further demonstrating the model's compatibility with existing experimental and modeling results on CaMKII activation (De Koninck and Schulman, 1998; Dupont et al., 2003; Chang et al., 2017). As the study of cellular information processing shifts from individual transduction pathways, toward the emergent properties of complex signaling networks, simple mathematical models are becoming indispensable tools for both experimentalist and theoreticians alike by providing a trade-off between detailed performance and a reduced description that facilitates system-level studies (Bornholdt, 2005; Kotaleski and Blackwell, 2010). Much in the way that the leaky-integrate and fire model has benefited the study of spiking neurons (Jolivet et al., 2004; Burkitt, 2006), the minimal switch model discussed in this paper will hopefully facilitate further study of complex kinase-phosphatase networks.

Results

A Bistable Switch Model

The following differential equation is an abstraction of a bistable molecular switch and was originally proposed as a model of genetic development by Lewis et al. (1977). This relatively simple model is a useful analytical tool to understand the general properties of bistable kinetic systems and captures the qualitative dynamics of more complicated models of CaMKII (Zhabotinsky, 2000) (Figure 1). Although the model interpretation and results presented here are centered on CaMKII and synaptic plasticity, the reader is encouraged to consider the broader implications for instantaneous frequency coding with other molecular switches, such as mitogen-activated protein kinases (Xiong and Ferrell, 2003; Thomas and Huganir, 2004).

$\begin{array}{l}\frac{dy}{dt}=k_{0}s-k_{1}y+\frac{k_2{y}^n}{{k_3}^n+y^n}\end{array}$

In this formulation, the level of activated CaMKII (y) is stimulated by the presence of Ca²⁺ bound to CaM, s, which will be studied as a function of time. For simplicity, it's assumed that pulses of Ca²⁺ are bound upon cell entry, which is reasonable since CaM is found in large concentrations surrounding Ca²⁺ channels and has a strong affinity for Ca²⁺ (Chin and Means, 2000). Switch deactivation is directly proportional to the active CaMKII concentration at a rate k₁, representing the activity of protein phosphatases. Finally, once activated, CaMKII has the ability to cooperatively bind Ca²⁺- CaM and autophosphorylate its own subunits, which motivates the nonlinear, positive feedback term captured by the Hill equation, where k₂ and k₃ are the association and dissociation constants, respectively. In addition to phosphorylation among the twelve subunits of a single CaMKII molecule, the ability to exchange active subunits between distinct CaMKII enzymes may connect this simple interpretation to a total, large pool of activated subunits distributed over multiple molecules (Stratton et al., 2014). Due to physiological constraints, y,s,k₀,k₁,k₂,k₃≥0. In the following, this specific equation will be referred to as the full kinetic model.

The full kinetic model of Lewis et al. has been previously applied to bistable genetic networks (Lewis et al., 1977; Smolen et al., 1998; Zheng et al., 2011), transcriptional regulation (Heltberg et al., 2016; Kang et al., 2017), mitogen-activated protein kinases (Xiong and Ferrell, 2003), and incorporated into a larger phenomenological model of presynaptic plasticity (Oswald et al., 2002). Although insightful for their specific systems, these studies retain a large numbers of parameters that clutter analysis and obscure the generality of the results. Therefore, it is desirable to reduce the number of parameters and facilitate the following analysis by performing routine nondimensionalization. Let y = x·k₃, $r=\frac{k_1{k}_3}{k_2}$, $s=\frac{k_2}{k_0}c$ and $t=\frac{k_3}{k_2}\tau$, which, when substituted into the original equation and simplifying gives the reduced but dynamically equivalent form:

$\begin{array}{rc}\frac{dx}{d\tau }=c-rx+\frac{x^n}{1+x^n}& (1)\end{array}$

This article is interested in a time varying c ≡ c₀ + c_l(τ), where c₀ reflects residual cytosolic Ca²⁺, whose slow dynamics are treated as fixed on the fast timescales over which the local Ca²⁺ signal c_l(τ) fluctuates (Regehr, 2012). A timescale factor $T=\frac{k_3}{k_2}$, the quotient of the switch deactivation and activation parameters, will be reintroduced later in order to connect the switch dynamics to time in seconds and stimulation frequency in Hz. The parameter r represents the kinetics of CaMKII subunit dephosphorylation by protein phosphatases and scales with the factor T. Finally, for highly cooperative reactions, n = 2 is a reasonable approximation of the Hill function exponent (Edelstein-Keshet, 2005) and a convention maintained by all of the studies listed above. The following bifurcation analysis is illustrated for n = 2, which allows for exact analytical solutions (Figure 1 and Methods); however, the main results are then generalized to arbitrary n ∈ ℝ⁺, which is much more realistic and has important consequences for frequency coding.

Stability and Bifurcation Analysis

Although interested in frequency-driven fluctuations over short timescales (Box 1), we must first examine the bistable, long timescale equilibrium behavior of the model that defines the switch's low and high activation states (Equation 1; Figure 1). An important reason for reducing the number of model parameters above is to simplify the analysis of all the possible system behaviors as a function of only a few parameter values. Having selected n = 2, we now only need to consider the effect of varying r and c; depending on their values, we may have one, two or three equilibrium points (x*), where the rate of change of the switch $f(x)=c-rx+\frac{x^2}{1+x^2}$ is equal to zero. For example, consider the values r = 0.52 and c = 0.04 that support bistability: there are three fixed points, two of which are stable, as illustrated by the switch's potential function U(x) = −∫f(x)dx (Figure 1A). As r and c change, saddle node bifurcations can occur, resulting in the presence of only the high or low activation state. The corresponding bifurcation diagrams are displayed in Figure 1B; their derivation is found in the Methods section.

A key feature of bistability is the hysteresis effect, where the same value of a parameter may evoke different states depending on the history of activity. For example, as the Ca²⁺ signal c increases, x* grows larger until crossing the rightmost c_c, where a saddle node bifurcation occurs and the switch jumps up to the high activation state, as the low state disappears (Figure 1Bi). Now, as c decreases back into the bistable range, the high activation state is preserved, and only lost when c crosses below the leftmost value of c_c. This history dependent behavior is presumably central to sustained CaMKII activity on the order of seconds (Wang and Maler, 1998) (Box 1). A similar phenomenon occurs for the negative regulation parameter r (Figure 1Bii). The values of r_c and c_c are plotted parametrically as a function of the active switch in the bifurcation curves (Figure 1Biii). The bifurcation surface summarizes this information completely (Figure 1Biv).

Existence of Solutions Around Stable Equilibria

To date, studies of Lewis et al.'s full kinetic model have been restricted to static input and periodic forcing. It is of principal interest to characterize the model behavior in response to aperiodic forcing, in order to gain a more general, physiologically realistic understanding of frequency coding with molecular switches. In addition to potentially encoding frequency information into stable levels of activated switch for many seconds presynaptically (Wang and Maler, 1998), or minutes postsynaptically (Lisman et al., 2012), what about frequency coding on the order of milliseconds to seconds, which is associated with brief sequences of action potential-evoked Ca²⁺ inputs? In a neighborhood surrounding a stable activation state (a sub-state region), is there a unique solution for a given time varying input signal? This question is not trivial, since small changes in the initial conditions of a nonlinear system (i.e., past switch activity) may generate drastically different behavior. Understanding the relationship that determines whether solutions converge or diverge around a given steady state could provide valuable insight into the properties of bistable molecular switches.

In the following section, we now reintroduce the scale factor T, since we are interested in studying frequency in Hz and time (t) in seconds. As such, Equation 1 becomes

$\begin{array}{rc}T\frac{dx}{dt}=c(t)-rx+\frac{x^n}{1+x^n}& (2)\end{array}$

First, to establish the existence of solutions around the high and low switch states, consider Equation 2 and note that f explicitly depends on the time-varying forcing term, c(t) ≡ c₀ + c_l(t). The phosphatase activity r that can counteract the switch phosphorylation is treated as fixed. The function f(t, x(t)) is assumed to be Lipschitz continuous and well-defined within intervals of state space, y₋ ≤ x(t) ≤ y₊ satisfying the conditions f(t, y₋) >0 and f(t, y₊) < 0 for all t ∈ ℝ⁺ (recall Figure 1A), which traps solutions within these boundaries. For any given point in time, there exist boundaries (y₋, y₊) determined by the parameters r, c₀, and the input c_l(t); we refer to values of the activated switch falling within these trapping regions as sub-state solutions, that is, fast timescale fluctuations that occur around either the high or low stable activation states (Lisman et al., 2012).

For (c, r) corresponding to the bistable region of parameter space (Figure 1Biii), there exist two intervals, x(t) ∈ (y_{l −}, y_{l +}) and x(t) ∈ (y_{h −}, y_{h +}), each surrounding one of the stable equilibrium points (x*). Now, we wish to locate values for the low state (y_l−, y_l+) and high state (y_h−, y_h+), where the existence of local time-varying solutions can be established. This problem is intimately linked to bifurcation, since y_l+ and y_h−depend on the values of c and r. The choice of a lower bound for the interval that exists around the low activation state is y_l− = 0, since the physiological restriction c(t) ≥ 0 implies f(t, 0) > 0 for all t ∈ℝ⁺, ignoring the boring degenerate case of c(t) = x(t) = 0. The upper bound of the lower strip, y_l+, can be chosen as a value $x_u^{{}_{*}}-\Delta x$, left of the unstable equilibrium $x_u^{{}_{*}}$ where $f(t,x_u^{{}_{*}})=0$, such that $\Delta c+f(t,x_u^{{}_{*}}-\Delta x)<0$; this condition ensures that the system is not trivially displaced into the high activation state by a single Ca²⁺ pulse with amplitude Δc. For the high concentration strip (y_h−, y_h+), the lower bound y_h− is chosen as a value of x infinitesimally greater than $x_u^{{}_{*}}$, that is, $y_{h-}=x_u^{{}_{*}}+\epsilon$ for ε → 0. Since we have restricted r and c>c_c (leftmost; Figure 1B) to the bistable range, we know that f(t, y_h−) > 0. For the upper bound of the high activation strip, it is enough to note that for $x>x_h^{{}_{*}}$, f(t, x(t)) < 0 and, since we wish to maximize the width of the strip, we take x arbitrarily large, denoting this value by y_h+ = x_∞. During stimulation, if (c, r) drifts out of the bistable region of parameter space, a saddle node bifurcation occurs and only one interval exists; in this case, the bounds simply span the state space, y₋ = 0 and y₊ = x_∞.

By invoking the Cauchy-Peano theorem, we guarantee the existence of at least one sub-state solution for every initial condition found within the interval regions defined above, since the conditions on the sign of the derivative f(t, x(t)) define trapping regions. However, this theorem says nothing about

$\begin{array}{l}\frac{d}{dt}z(t)=\underset{h\to 0}{\lim }\frac{z(t+h)-z(t)}{h}\\ =\underset{h\to 0}{\lim }\frac{\left|u(t+h)-x(t+h)\right|-\left|u(t)-x(t)\right|}{h}\\ \le \underset{h\to 0}{\lim }\frac{\left|(u(t+h)-x(t+h))-(u(t)-x(t))\right|}{h}\\ =\underset{h\to 0}{\lim }\frac{\left|(u(t+h)-u(t))-(x(t+h)-x(t))\right|}{h}\\ =\text{sgn}[u(t)-x(t)]·\frac{d}{dt}(u(t)-x(t))\\ =T^{-1}\text{sgn}[u(t)-x(t)]\left(c(t)-r·u(t)+\frac{u^n(t)}{1+u^n(t)}-\left(c(t)-r·x(t)+\frac{x^n(t)}{1+x^n(t)}\right)\right)\\ =T^{-1}\text{sgn}[u(t)-x(t)]·\left(-r·(u(t)-x(t))+\frac{u^n(t)-x^n(t)}{(1+u^n(t))(1+x^n(t))}\right)\\ =T^{-1}\text{sgn}[u(t)-x(t)]·(u(t)-x(t))·\left(-r+\frac{u^n(t)-x^n(t)}{(u(t)-x(t))(1+u^n(t))(1+x^n(t))}\right)\\ =T^{-1}\left|u(t)-x(t)\right|·\left(-r+\frac{{\displaystyle \sum _{i=1}^n{u}^{i-1}(t)·x^{n-i}(t)}}{(1+u^n(t))(1+x^n(t))}\right)\text{for}n\in {\text{Z}}^{+}\\ =T^{-1}z(t)·\left(-r+p(u,x,n)\right)\end{array}$

whether solutions starting at different initial conditions will converge to a unique, stimulus-driven response that tracks changes in the Ca²⁺ signal.

Uniqueness of Sub-State Solutions

As motivation for the following results, Figure 2A shows an example switch response to an 8 Hz Poisson pulse sequence, which is convolved with an alpha function filter (30 ms, Methods), then normalized to the signal's maximum and scaled by Δc = 0.5 to create an example input signal, which the switch tracks closely. Note, in this simulation, the alpha-function kernel was specifically chosen to be 30 ms based on literature values for the time course of local synaptic Ca²⁺ signals (Sinha et al., 1997; Sabatini et al., 2002; Graupner and Brunel, 2012). Due to our interest in the fast timescales associated with short sequences of input pulses (100s of milliseconds), we assume that the Ca²⁺ pulse size is fixed on this timescale, which is a reasonable approximation for hippocampal spiking frequencies less than 15 Hz (Tank et al., 1995). This distinction between short and long timescales provides a hypothetical means for the system to be less sensitive to variations in the Ca²⁺ pulse size and the resulting frequency-intensity coding ambiguity [(Zhabotinsky, 2000; Pinto et al., 2012); see Box 1]. This could allow for more accurate representations of instantaneous frequency over short time periods, compared to long timescale frequency coding where input history, as well as additional adaptive and homeostatic processes may substantially adjust Ca²⁺ signaling.

FIGURE 2

Figure 2. Switch activity fluctuates with instantaneous input frequency. (A) Motivating example: switch response to an 8 Hz Poisson sequence of input pulses, convolved with an alpha function kernel to create a signal, c(t) = c₀+c_l(t). The switch's fluctuations track changes in the input frequency (n = 1.6, r = 0.61, c₀ = 0.04, and T = 0.01). (B) The example function p(u, x, 2) from the uniqueness proof achieves a maximum of 0.65; r must exceed this value to guarantee absolute convergence of the switch to a unique frequency-driven solution. (C) Initial conditions: u(0) = 1.7 and x(0) = 0.1. The value of T affects time-to-convergence between solutions and frequency filtering. From empirical studies, T ≤ 0.01 (Coultrap and Bayer, 2012). (D) (i) For r = 0.54 < r_c, sufficiently high frequency Ca²⁺ pulses (bursts) cause transitions from the low to high state (illustrated for n = 2). By adjusting c₀ to take advantage of hysteresis, the cell can control whether or not it is sensitive to these burst-induced up states. The first two pulses (< 10 Hz), where c₀ = 0, do not result in hysteresis, whereas the high frequency 10 Hz inter-pulse interval (right panel black arrow) with c₀ augmented to 0.04 does; note that neither static value can generate the upstate alone without sufficient input (e.g., 5 Hz, left panel black arrow). The switch response differs during the transition between low and high states, but once settled around a given state gives good agreement (gray shading; the two example curves are compared by choosing an offset of 0.92 that minimizes the Euclidean distance between them). (ii) Top Simulation for r = 0.59 and n = 2, where x(t) has c₀ = 0.02 and thus cannot support bistability, vs. u(t) with c₀ = 0.04, which, when driven by input, traps the solutions around the high activation state through hysteresis. Under these conditions, convergence cannot occur. Bottom The absolute value of the difference between the relative changes in u(t) and x(t) induced respectively by the common input frequencies (Δ, see Results for details) plotted as a function of r; the maximum discrepancy of 0.046 is found at r = 0.59 and represents a small fraction of the total activated switch concentration. (E) In general, the exponent n≠2 in real biological systems. Interestingly, n = 1.55 is a minimum for the maximum value of the class of functions p(u, x, n) in the uniqueness proof. This is remarkably close to the empirical best-fit value of 1.6 their ± SEM or SD reported by De Koninck and Schulman for presynaptic α CaMKII (De Koninck and Schulman, 1998).

We now establish the stability and uniqueness of solutions for distinct initial conditions within a given interval of state space. Consider a general interval (y₋, y₊), where x(t) is a solution to Equation 2 with initial condition x₀ ∈ (y₋, y₊). Assume there is another solution, u(t), with a different initial condition u₀ ∈ (y₋, y₊). Writing z(t) =|u(t)−x(t)| and first assuming n is a positive integer, we see that

The expression p(u, x, n) achieves maximal values at intermediate switch levels that separate the low and high states of activation. Now, consider p(u, x, n) for the special case of n = 2 used in the bifurcation analysis; in this case, $p(u,x,2)=\frac{u+x}{(1+u^2)(1+x^2)}$, which is plotted in Figure 2B. Setting the partial derivatives of the function to zero and solving for u and x, yields a critical point: $(u,x)=(\frac{\sqrt{3}}{3},\frac{\sqrt{3}}{3})$. Substituting this into p gives a global maximum of $\frac{3\sqrt{3}}{8}$. Since $\frac{d}{dt}z(t)\le T^{-1}(-r+\frac{3\sqrt{3}}{8})z(t)$ for all t, we can apply Grönwall's inequality, which gives us the following:

$z(t)\le e^{T^{-1}}\underset{0}{\overset{t}{{{\displaystyle \int }}^{\text{}}}}(-r+\frac{3\sqrt{3}}{8})ds$

Substituting the expression for z(t) and solving this integral exponent yields,

$|u(t)-x(t)|\le e^{-T^{-1}(r-\frac{3\sqrt{3}}{8})t}$

and, as t → ∞, we have

$\begin{array}{l}0\le \underset{t\to \infty }{\lim }|u(t)-x(t)|\le \underset{t\to \infty }{\lim }{e}^{-T^{-1}(r-\frac{3\sqrt{3}}{8})t}\hfill \end{array}$

For $r>\frac{3\sqrt{3}}{8}$ (≈ 0.65), we obtain

$\begin{array}{l}0\le {\displaystyle \underset{t\to \infty }{lim}}|u(t)-x(t)|\le 0\end{array}$

By the squeeze theorem we conclude that |u(t) − x(t)| → 0 as t → ∞. Therefore, a unique frequency-driven solution exists and is independent of the initial conditions within the bounded interval. The time taken to converge to the unique solution is inversely proportional to $T=\frac{k_3}{k_2}$ (Figure 2C). The parameter value T = 0.01 s was chosen here for our specific example switch, CaMKII, whose dissociation constant (k₃) has been experimentally determined to be at least 100-fold smaller than the activation constant (k₂) that governs the rate of autophosphorylation (Coultrap and Bayer, 2012). Unlike the larger value of T = 0.1 s, T = 0.01 permits quick convergence and reliable encoding for the action potential frequencies characteristic of hippocampal CA3-CA1 synaptic input (approximately 1–15 Hz) (Mizuseki et al., 2012). Smaller values of T permit rapid convergence and more accurate frequency coding, but may become overly sensitive to temporary lulls in activity when c briefly drops below the leftmost critical value c_c (recall Figure 1Bi).

It should be noted that r > 0.65 is an absolute guarantee of convergence to a unique frequency driven solution; but, from the bifurcation analysis (Figure 1Biii; Methods), we know that bistability does not exist for this value of r. However, in general, only −r + p(u, x, 2) < 0 is required, which, for low and high concentrations of activated switch, is obtained at smaller values of r that do support bistability. In fact, p(u, x, n) only exceeds the r value briefly during state transitions as it moves through the unstable equilibrium. Although a unique encoding of sub-state solutions can still exist for smaller r values around either the high or low state, convergence about the low activation state is now vulnerable to perturbation by short Ca²⁺ inter-pulse intervals, thus acting as a high frequency event (burst) detector through induction of high switch activation. For example, experiments show that high frequency hippocampal activity (>15 Hz) causes successive Ca²⁺ pulses to accumulate (Tank et al., 1995), which could effectively boost c₀ and serve to promote burst detection by transiently maintaining the high activation state via hysteresis (Figure 2Di). In theory, this dynamic burst threshold (the separatrix) is sensitive to recent levels of activation, and could be purposefully modulated by the cell through dynamic regulation of the parameters r and c₀ (Li et al., 2012). To restore the low state, the cell simply needs to adjust c₀ to fall below the leftmost critical value c_c. The bottom panel of Figure 2Di illustrates that fluctuations around the high- and the low-stable states still yield nice agreement in their response to a given input frequency. Of course, during the state transition itself, the switch response can differ largely but once solutions are settled around their respective stable states the model appears to give good agreement (gray shading; the two example curves are compared by choosing an offset of 0.92 that minimizes the Euclidean distance between them).

When bistability is supported, the model response cannot always converge to an absolute level of activated switch, as illustrated in the top panel of Figure 2Dii; however, the fluctuations about the distinct stable states appear to be similar, as in Figure 2Di. To examine this idea further, repeated simulations of the model were performed, where x(t) has an associated residual Ca²⁺ level of c₀ = 0.02 and thus does not support bistability, vs. u(t) with c₀ = 0.04, which can trap the solution around the high activation state through hysteresis (Figure 2Dii, top). As was the case in Figure 2C, the same random spike sequences are used for x(t) and u(t) on each trial. For each inter-pulse interval of the repeated simulations, the change in the level of activated switch was computed as the difference between the switch activity sampled at the time of an input pulse and the subsequent maximum switch response that occurred before the next pulse. For each successive, shared inter-pulse interval, these differences, Δx and Δu, were determined separately for x(t) and u(t), then subtracted from each other for each 100 s trial, containing an average of 797 pulse intervals (8 Hz Poisson process). This was repeated 10 times for each parameter set and the composite mean of the absolute value of the difference between the change in the two solutions, Δ = |Δu−Δx|, was determined as a function of r (Figure 2Dii, bottom). The maximum discrepancy between Δx and Δu, 0.046, occurs at r = 0.59 (used for Figure 2Dii, top) and is at least an order of magnitude less than typical values achieved in the low activation state. These results suggest that the relative change in switch activation about a stable state is generally quite consistent.

Realistically, the Hill function exponent n need not be restricted to integer values, which is unlikely in real biological systems. Thus, in the above proof, the expression p(u, x, n) is now left as $\frac{u^n-x^n}{(u-x)(1+u^n)(1+x^n)}$ for n ∈ ℝ⁺, since there is no longer a closed form expression for the factorization of the numerator by the term u − x. The function p(u, x, n) has critical points at u = x, which occur at an apparent discontinuity due to the factor u − x in the denominator. However, assessing the limit as the difference between x and u becomes infinitesimally small, making the change of variable u = x + h as h → 0, and recognizing the limit definition of the power rule for differentiation, yields an expression for the maximum of p(u, x, n) for all u, x, n ∈ ℝ⁺:

$\begin{array}{l}\max [p(u,x,n)]=\underset{u\to x}{\lim }p(u,x,n)\\ =\underset{u\to x}{\lim }\left[\frac{u^n-x^n}{(u-x)(1+u^n)(1+x^n)}\right]\\ =\underset{h\to 0}{\lim }\left[\frac{{(x+h)}^n-x^n}{(x+h-x)}\right]\\ ·\underset{h\to 0}{\lim }\frac{1}{(1+{(x+h)}^n)(1+x^n)}\\ =\frac{d}{dx}\left(x^n\right)·\frac{1}{{(1+x^n)}^2}\\ =\frac{n{x}^{n-1}}{{(1+x^n)}^2}\end{array}$

For each value of the exponent n, the global maximum of this expression is determined for all x ∈ ℝ⁺, and plotted (Figure 2E). Ignoring the highly uncooperative reaction exponents of n < 0.012, the global minimum of the class of functions p(u, x, n) is found at n = 1.55. Fascinatingly, the empirical αCaMKII data reported by De Koninck and Schulman was fit by a Hill function with an exponent of 1.6 (De Koninck and Schulman, 1998; Dupont et al., 2003). This intriguing match between their experiment and the model's theory suggests that αCaMKII's activation function may operate with this particular exponent as it provides the minimum level of negative regulation r required to maintain absolute convergence of unique input driven switch activity in the low activation state, even for intermediate levels of the switch response occurring just left of the unstable equilibrium (Figure 1A), where r must be much stronger to guarantee unique solutions (Figure 2B). As we will see in the following section, the value of n = 1.6 has additional benefits for amplifying the frequency response of weak calcium fluctuations in the presence of noise.

Molecular Switches and Stochastic Resonance

If Equation 2 is to capture actual molecular switch behavior in vivo, then we must understand frequency coding in the presence of biological noise. Given our interest in synaptic information transfer, it is natural to ask whether noise can improve the switch's frequency coding ability through stochastic resonance and how different combinations of our main parameters (for example the value of n) could potentially affect this phenomena. In particular, does the value n = 1.6 confer benefits for frequency coding? The results presented in this section are generated by Equation 2 with additive Ornstein-Uhlenbeck noise, η(t), which evolves according to the stochastic differential equation

$\begin{array}{l}\frac{d\eta }{dt}=-\frac{\eta }{{\tau }_{\eta }}+\xi (t)\end{array}$

where ξ(t) is bounded Gaussian noise, N(0, 1), whose amplitude is scaled by a parameter σ. The simplest interpretation is that there is some weak noise in the Ca²⁺ signal amplitude, which might arise from stochastic channel dynamics. The choice of the time constant τ_n is based on previous studies of noisy microdomain Ca²⁺ fluctuations, where an upper bound for the autocorrelation time was determined to be approximately 10 ms (von Wegner et al., 2014; Weinberg and Smith, 2014). This choice has the added benefit of matching our switch time constant T, should we instead assume the noise is inherent to switch activation, as well as matching a typical value for the membrane time constant of spiking neurons, whose noisy membrane potential fluctuations might influence the activity timescales of voltage-gated Ca²⁺ channels.

Figure 3A shows the power spectrum (P_c) of a weak sinusoidal Ca²⁺ oscillation, c = c₀+αsin(2πφt), where c₀ = 0.04, α = 0.02 and φ = 2 Hz, which was selected based on the mean action potential frequency associated with the CA3 and CA1 regions of the hippocampus (Csicsvari et al., 2000). As expected, the noisy switch oscillates at the frequency φ, reflected in its power spectrum (P_x). Very recently, the full kinetic model of Lewis et al., studied under the context of genetic regulation with n = 2, has been shown to produce the stochastic resonance effect (Kang et al., 2017), which is confirmed here for the dimensionally reduced model (Equation 2; Figure 3B). As σ increases from 0, frequency transfer, measured as the ratio of the switch power to signal power at φ, dips slightly and then improves dramatically, achieving a maximum at 0.29, followed by a quick decrease as the noise becomes dominant. When changing the exponent from n = 2 to n = 1.6, this spectral amplification becomes significantly larger, further suggesting that presynaptic α CaMKII functions as an important frequency decoder and that the exponent n = 1.6 may have evolved to fulfill this purpose. The reader should note that, for fair comparisons sake, r = 0.65 and r = 0.61 were selected respectively for n = 2 and n = 1.6 based on values obtained from Figure 2E, but this effect is qualitatively robust to changes in r and φ. Setting n = 1.6 also shifts the optimal noise strength to a substantially lower value, 0.09, which has the putative benefit of harnessing stochastic resonance and enhanced frequency representations for low intensity Ca²⁺ signal noise.

FIGURE 3

Figure 3. Frequency coding with noisy switches. (A) The switch model driven by a weak sinusoidal signal, c(t) = c₀+αsin(2πφt), with c₀ = 0.04, α = 0.02, φ = 2 Hz, and additive noise, η(t), whose intensity is scaled by the parameter σ and evolves according to τ_η = 0.01. The switch amplifies the frequency content of the input, as shown by its power spectrum P_x relative to the signal's, P_c. (B) Top: For n = 2, the ratio of switch power to signal power at φ is plotted as a function of the noise intensity σ, achieving a maximum at 0.29, that is, the switch displays stochastic resonance (SR). The value of σ that promotes optimal frequency transfer is denoted by σ_p. Bottom: For n = 1.6, there is substantially larger gain in the SR effect, and σ_p shifts to 0.09. (C) σ_p is plotted as a function of (n, r) and (c₀, r), illustrating the presence or absence of SR. (D) For n = 1.6, stochastic switch simulations produce bimodal (e.g., σ = 0.01) or unimodal (e.g., σ = 0.035) activation around the low state (left column; r = 0.61, c₀ = 0.04) and the high state (right column; r = 0.54, c₀ = 0.04) (Box 1). Within each sub-state region, the input is uniquely encoded. (E) As model validation, the pulse duration (ms) and frequency experiments of De Koninck and Schulman were simulated (n = 1.6, r = 0.61, and T = 0.4), qualitatively capturing their results, as well as the results of a follow-up model (Dupont et al., 2003). The reader should note the ambiguity in autonomously activated (long timescale) switch activity, based on input duration and frequency (Pinto et al., 2012). (F) To validate the model against short timescale CaMKII data, a 1 Hz pulse train generates a calcium signal that drives the switch model: the top plot shows the entire period of stimulation and the decay of the switch after cessation of the stimulus (n = 1.6 and T = 0.1). These simulations qualitatively capture the exciting new experimental data of Chang et al. (2017). Note that, instead of a fixed value of negative regulation, r now linearly increases from 0.58 to 0.67 over the course of the stimulus, providing a potential explanation of the slow decay in the plateau switch concentration seen in the data but not in the model for a fixed value of r. The bottom plot zooms in on the first 15 s and also plots the model prediction for a T286A mutant form of CaMKII, which prevents the nonlinear switch activation. Naturally, this results in a much weaker response compared to the wild-type (WT), as seen in the data (Chang et al., 2017).

The model results of Kang et al. (2017) depend on a full complement of parameters, which begs the question of whether stochastic resonance is a generic feature of the model switch or whether the effect is only significant for a certain range of the parameters. The dimensional reduction of the switch model performed here allows this question to be easily addressed as a function of the parameters c₀, r and n. Figure 3C shows that the parameter r has significant influence over the value of σ that produces optimal spectral amplification and that, for some combinations of c₀, r and n, the stochastic resonance effect disappears completely. The presence or absence of stochastic resonance may prove useful for deducing parameter ranges of molecular switches in vitro and in vivo. Furthermore, these noise fluctuations may generate unimodal (e.g., σ = 0.035) or bimodal (e.g., σ = 0.01) distributions of switch activation (Figure 3D, n = 1.6), which provides another experimentally testable prediction for αCaMKII, given that the switch state could control neurotransmitter release (see Box 1) and thus explain multimodal distributions of excitatory postsynaptic potential amplitudes (Larkman et al., 1992). The occupation of the low state (Figure 3D, left) vs. the high state (Figure 3D, right) depends on the level of negative regulation r and whether c₀ can support hysteresis: the parameter choices for the left column of Figure 3D do not support bistability (r = 0.61, c₀ = 0.04) and the switch fluctuates around the low activation state. The right column of Figure 3D does support bistability (r = 0.54, c₀ = 0.04) and input activity quickly drives high switch activation, while hysteresis ensures the switch stays within this state. Stochastic simulations for Figure 3D were performed by including additive Ornstein-Uhlenbeck noise, as described above. Further detail can be found in the Methods section.

Bridging Short Term Dynamics With Long Timescale Switch Activation

A potential caveat of the bistable switch model is that, even in the high activation state, the population of phosphorylated units (x) are still subject to the phosphatase activity (r). Equation 2 places difficult constraints on cells for long-timescale activation: if c₀ and r are not controlled carefully, the high activation state can be lost. Although high activation levels may only be short lived in vivo, it is important to establish a potential connection between the current model and existing theories of long timescale activation (Box 1; Introduction). Equation 2 effectively represents all of the phosphorylated subunits in a population of CaMKII molecules (each having twelve phosphorylation sites). When one of these dodecamers becomes fully phosphorylated, it could effectively become impervious to negative regulation by the phosphatases, since any cleaved subunit could immediately be re-phosphorylated by its neighboring subunits and the enzyme can be shielded by its interactions with downstream targets (e.g., an NMDA receptor subunit) (Lisman et al., 2012; Urakubo et al., 2014). Until now, the work presented here has ignored this potentially important feature of CaMKII, since the actual biological relevance of autonomous activation is still in question (Box 1). Therefore, to connect the short term dynamics to long timescales, we introduce a new variable (X) to represent the level of autonomously activated switch that might persist after the stimulus has been removed, even when Ca²⁺ levels drop below the leftmost critical value c_c that supports hysteresis (Figure 1Bi). X is calculated from Equation 2 by using Equation 3, explained below.

Motivated by the work of Pinto et al. (2012) (Box 1), let us assume that the total amount of autonomously activated switch (X) is simply proportional to the average amount of Ca²⁺ input, which is determined by pulse amplitude, duration and frequency. As seen in Figure 2, this value is reflected by the amount of activated switch x(t) over the duration of the stimulus, Δt. Therefore, let X be the temporal average of x(t)

$\begin{array}{cc}X=\langle \omega x(t)\rangle =\frac{\omega }{\Delta t}\underset{0}{\overset{\Delta t}{{{\displaystyle \int }}^{\text{}}}}x(t)dt& (3)\end{array}$

The biological interpretation is as follows: at a given moment in time there is some likelihood for individual dodecamers to transition to the fully autonomous, phosphorylated switch state or bind to a downstream target. These autonomous elements accumulate over time. For simplicity, a fixed basal rate of transition of a given molecule to the fully autonomous state, ω, is assumed.

As a validation of the model's ability to produce CaMKII-like behavior over long timescales, the essence of De Koninck and Schulman's experimental results (De Koninck and Schulman, 1998) and the model of (Dupont et al., 2003) are both captured qualitatively by Equations 2 and 3 (Figure 3E). Note that this result was generated using Equation 2 and 3, but does not include Ornstein-Uhlenbeck noise given the synthetic and controlled nature of the original experiment (De Koninck and Schulman, 1998). The timescale factor T was set on the order of 10⁻¹ s, which may reflect altered kinetics under the artificial conditions of the experiment, or the need for further refinement of the model presented here. For instance, the proportion ω is expected to grow larger as more of the subunit population becomes phosphorylated and cooperative activation grows stronger (Meyer et al., 1992; Chao et al., 2010), leading to an increased likelihood for individual dodecamers to transition to the fully autonomous state. This is expected to improve the reproduction of De Koninck and Schulman's results by flattening the curves at lower frequencies and steepening them at higher frequencies (De Koninck and Schulman, 1998). Future work should seek to determine ω(x), with the hopes of identifying reduced representations of strongly nonlinear CaMKII activation. In general, ω could also depend on the interaction of activated CaMKII with downstream targets; for example, in vitro evidence suggests constitutive CaMKII activation (hysteresis) requires interaction with an NMDA receptor peptide (Urakubo et al., 2014).

Given our interest in the short timescale behavior of the CaMKII switch, Equation 2 was further validated against recent experiments that used fluorescent life-time imaging microscopy to measure CaMKII activity with millisecond precision in neurons responding to glutamate uncaging (Chang et al., 2017). Figure 3F shows 1 Hz stimulation (30 pulses) that drives a dynamic calcium concentration, modeled as a first-order exponential decay with a time constant of 200 milliseconds. This signal (c; not shown) is used to drive the switch model with exponent n = 1.6. The value of T was set to 0.1 for this simulation; as an aside, T may be impacted by the mean frequencies experienced at the synapse over long timescales and can be modified according to the specifics of the system (e.g., pre- vs. post-synaptic), as well as further influenced by experimental conditions. These factors might explain the difference used to fit the Chang et al. (2017) from the value of T = 0.01 inferred from Coultrap and Bayer (2012). To account for the slow decay of plateaued switch activity (Chang et al., 2017), r was made to linearly increase for a small range over the course of stimulation. The justification for this is as follows: the protein phosphatase calcineurin has a much higher affinity for Ca²⁺ ions than calmodulin, and will slowly strip Ca²⁺ away from calmodulin that has yet to bind to CaMKII. Thus, it seems possible that the accumulation of activated calcineurin over the course of prolonged stimulation may exert a growing, adaptive effect on the switch. Although the linear change in r was chosen for simplicity, it's likely more complicated and may account for minor discrepancies between the model and the actual data of Chang et al. (2017). Finally, despite the fact that constitutively active CaMKII is unlikely to occur in vivo, the bistable-associated nonlinearity that governs its dynamics is still essential to its function. By repeating their experiment with a mutated form of CaMKII that cannot autophosphorylate (T286A), Chang et al. provide evidence that T286 phosphorylation is essential for the optimal integration of Ca²⁺ signals by boosting the activation of the switch and slowing its decay. The bistable model nonlinearity in Equation 2 reflects this cooperative autophosphorylation – removing it results in a much weaker response (Figure 3F), as seen in the data (Chang et al., 2017). In addition to enabling the induction of long term plasticity postsynaptically, it is hypothesized that this nonlinearity is also essential for the switch to sequester cytosolic calcium and thus regulate neurotransmitter release at the presynaptic terminal.

Discussion

A main goal of this study was to extend the frequency coding idea of De Koninck and Schulman (1998) in a generic switch model that captures the qualitative behavior of CaMKII, but focuses on fast timescale dynamics instead of slow timescales (Box 1). The model presented here may help to reconcile contradictory perspectives of CaMKII function (De Koninck and Schulman, 1998; Pinto et al., 2012) and suggests dual streams of information transfer that are temporally multiplexed: over short timescales, where the size and duration of the Ca²⁺ pulse are more stable (Tank et al., 1995), the molecular switch can act as an encoder of instantaneous frequency information (e.g. Figure 2A) and function to bidirectionally regulate transmitter release at synapses through a combination of enzymatic and non-enzymatic activity (summarized in Box 1). Over longer timescales, the model switch integrates overall signal intensity, which could dictate long term changes in synaptic strength and is dependent on multiple factors such as slow Ca²⁺-induced Ca²⁺ release (affecting c₀) (Sharma and Vijayaraghavan, 2003; Catterall and Few, 2008), the size of the Ca²⁺ pulse, its duration and the mean frequency of stimulation (Figure 3E).

Although the present work is a very preliminary investigation of the role of molecular switches in the processing of information in the brain, it provides some testable predictions for synaptic physiologists: establishing the presence of both bimodal and unimodal synaptic release that depends on αCaMKII and noise, as well as characterizing the hypothesized real-time modulation of release probability at central synapses by αCaMKII in response to natural, aperiodic stimulation patterns (specifically detection of bursting events). Finally, of particular interest, is the putative role of αCaMKII in the regulation of synchronous discharge probability and duration, as well as the propagation of CA3 oscillations into the CA1 area (Hinds et al., 2003). A more complete study including a coupling of the subcellular switch dynamics with those of the neuron and synapse will be published in the future.

Fascinatingly, the Hill function exponent of approximately 1.6 is not unique to CaMKII; calcium sensors within hippocampal basket cells display cooperative binding that also reflects the value 1.6, inferred from measurements of post-synaptic currents (Debanne et al., 2013). The mitogen-activated protein kinase (MEK1) is reported to have a Hill function exponent of approximately 1.7 (Ferrell and Ha, 2014), further hinting at the generality of the switch model. In general, bistable molecular switches such as CaMKII, are a conserved feature of cell signaling networks and generate combinatorial power in their collective action (Ferrell, 1997; Bhalla and Iyengar, 1999; Brandman et al., 2005). As previously described, stacking kinase pathways leads to an increase in the effective cooperative binding (described by the Hill function exponent) (Ferrell, 1997); for example, the extracellular signal-regulated kinase 2 (ERK2), which lies downstream of MEK1 is reported to have an approximately three-fold larger Hill coefficient of 4.9 (Ferrell and Ha, 2014), which is also associated with the famous Calyx of Held synapse (Debanne et al., 2013).

Due to their complex kinetics and network interactions, switch models are typically formulated by parameterized systems of differential equations that are not ideal for deeper mathematical analysis. It is proposed that the simple model described by Equation 2 can capture the core essence of molecular switches, much in the way that the leaky-integrate and fire model has been a successful abstraction of neuronal spiking activity, providing a trade-off between performance and a reduced description that facilitates network studies (Burkitt, 2006; Jolivet et al., 2008). This idea is supported by the inclusion of Equation 2 in an existing phenomenological model of feedback-driven synaptic plasticity, using the conventional exponent of n = 2 (Oswald et al., 2002). The relative simplicity of the switch model and its application to diverse signaling pathways make it a useful framework for further theoretical and experimental investigations into signaling networks, synaptic plasticity and cellular computation.

Methods

Bifurcation Analysis

The first step of the bifurcation analysis is to find the equilibrium points. Setting n = 2, we rewrite Equation 1 as,

$\begin{array}{l}\frac{dx}{dt}=g(x)-h(x)\end{array}$

where $g(x)=\frac{x^2}{1+x^2}$ and h(x) = rx − c. The fixed points occur when g(x) − h(x) = 0, which amounts to finding the solutions of the polynomial −rx³ + (c + 1)x² − rx + c = 0. First, fix c and examine the effects of varying r. When c = 0, x = 0 is a fixed point, and, for a particular range of r, there exists two other positive valued fixed points, given by the roots of −rx² + x − r = 0. The critical value of the parameter r, denoted by r_c is found by setting g(x) = h(x) and g′(x) = h′(x), which, when solved, gives $r_c=\frac{x}{1+x^2}=\frac{2x}{{(1+x^2)}^2}$. Three values of x satisfy this equality: −1, 0, and 1. Since we are not considering negative values of x, we have two critical points, r_c = 0 and $r_c=\frac{1}{2}$. Therefore, when c = 0, the system is bistable for $0<r<\frac{1}{2}$. For c>0, r can be larger than $\frac{1}{2}$ while still preserving bistability (as in Figure 1A). We know r_c occurs when h(x) = g(x) and h′(x) = g′(x); therefore, when h(x)>g(x) we lose a fixed point through a saddle node bifurcation. Forx>0, the maximum of g(x) is found at $x=\sqrt{\frac{1}{3}}$ which gives ${max}_x\left[g(x)\right]=\frac{3\sqrt{3}}{8}$. Therefore, when $r>r_c=\frac{3\sqrt{3}}{8}$, only one fixed point exists.

Now, we are interested in fixing r and examining the effects of varying c. To find c_c we set g(x) = h(x) and g′(x) = h′(x), which gives $r=\frac{2x}{{(1+x^2)}^2}$ and $c_c=rx-\frac{x^2}{1+x^2}$. Substituting the first expression into the second, we get $c_c=\frac{x^2(1-x^2)}{{(1+x^2)}^2}$. We differentiate with respect to x in order to locate the maximum value for c_c; $0=\frac{2x(1-3x^2)}{{(1+x^2)}^3}$. This gives x = 0 and $x=\sqrt{\frac{1}{3}}$, which corresponds to c_c = 0 and $c_c=\frac{1}{8}$. When c>c_c, only one fixed point exists for all values of r. For a fixed value of r that supports bistability, as c increases from 0 and crosses a critical value (c_c), the fixed point x* will jump up to the high amplitude branch. If c is now decreased, the fixed point remains on the high amplitude branch even as c becomes smaller than the corresponding c_c. This hysteresis effect permits switch activation to remain as the transient Ca²⁺ signal subsides, consistent with the findings from synaptic plasticity experiments (Box 1). Using the expressions derived for the critical values of r_c and c_c, we plot them parametrically as functions of x (Figure 1Biii). Saddle node bifurcations occur all along the boundary of these curves, it is here we find the values of r and c for which only two fixed points occur. Crossing each branch results in a pairwise collision and disappearance of two fixed points. Note where the bifurcation curves meet tangentially, $(c,r)\to (\frac{1}{8},\frac{3\sqrt{3}}{8})$, here we observe a co-dimension two bifurcation; beyond this point there is only one fixed point and the distinction between low and high activation states is blurred (Figure 1Biii).

Computational Specifications and Miscellaneous Details

Simulations were solved using the 4th order Runge-Kutta method, with the exception of the Ornstein-Uhlenbeck noise, which was solved using the stochastic Euler method (time step of 1 ms in all cases). All simulations were performed using custom code, available upon request to the author, and were implemented on a Linux machine running Ubuntu 16.04 with an Intel core i7-6700 CPU, 3.4 GHz processing speed, and 62 GB of RAM.

Pulse train sequences {t_i} were convolved with the filter $t·e^{-(t-t_i)/{\tau }_c}$, whose decay constant τ_c was set to 30 ms, reflecting an accommodation of both pre- and post-synaptic calcium decay values from the literature that range from 15 to 43 ms (Sinha et al., 1997; Sabatini et al., 2002; Graupner and Brunel, 2012). The resulting input signal was normalized to the maximum value and then scaled by Δc. The decay value is closely related to the input frequencies typical of a given synapse and the definition of what constitutes a high frequency event in the system, since for events occurring faster than the decay, Ca²⁺ accumulates quickly, driving the switch into the upstate. The putative burst detector will work for different τ_c, but may require a different set of corresponding switch parameters, range of stimulation frequencies and pulse amplitudes.

Histogram bin sizes for Figure 3D were set using the Freedman-Diaconis method (Freedman and Diaconis, 1981).

Author's Note

Bistable molecular switches can decode cellular inputs into distinct high- or low-states of persistent enzymatic activity. Although this on-off, “digital” perspective is valuable for long timescales, I suggest that short timescale fluctuations of switch activity around either stable state acts as an analog signal that reliably encodes instantaneous input frequency. A minimal model and theory make predictions about the molecular switch CaMKII, synaptic plasticity and burst detection.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

I would like to thank Lin Wang for useful discussions and guidance during the early stages of the project, as well as Leonard Maler, Richard Naud, and Jean-Claude Béique for proofreading the manuscript and providing valuable comments. This work was funded by training grants awarded to SC from both the National Sciences and Engineering Research Council of Canada and the Canadian Institute of Health Research. A preprint of the manuscript is available on bioRxiv (Clarke, 2018).

References

Berridge, M. J., Bootman, M. D., and Roderick, H. L. (2003). Calcium signalling: dynamics, homeostasis and remodelling. Nat. Rev. Mol. Cell Biol. 4, 517–529. doi: 10.1038/nrm1155

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhalla, U. S., and Iyengar, R. (1999). Emergent properties of networks of biological signaling pathways. Science 283, 381–387.

PubMed Abstract | Google Scholar

Bornholdt, S. (2005). Less is more in modeling large genetic networks. Science 310, 449–451. doi: 10.1126/science.1119959

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandman, O., Ferrell, J. E., Li, R., and Meyer, T. (2005). Interlinked fast and slow positive feedback loops drive reliable cell decisions. Science 310, 496–498. doi: 10.1126/science.1113834

PubMed Abstract | CrossRef Full Text | Google Scholar

Burkitt, A. N. (2006). A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input. Biol. Cybernetics 95, 1–19. doi: 10.1007/s00422-006-0068-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Catterall, W. A., and Few, A. P. (2008). Calcium channel regulation and presynaptic plasticity. Neuron 59, 882–901. doi: 10.1016/j.neuron.2008.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, J. Y., Parra-Bueno, P., Laviv, T., Szatmari, E. M., Lee, S. R., and Yasuda, R. (2017). CaMKII autophosphorylation is necessary for optimal integration of Ca(2+) signals during LTP induction, but not maintenance. Neuron 94, 800–808 e4. doi: 10.1016/j.neuron.2017.04.041

CrossRef Full Text | Google Scholar

Chao, L. H., Pellicena, P., Deindl, S., Barclay, L. A., Schulman, H., and Kuriyan, J. (2010). Intersubunit capture of regulatory segments is a component of cooperative CaMKII activation. Nat. Struct. Mol. Biol. 17, 264–272. doi: 10.1038/nsmb.1751

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapman, P. F., Frenguelli, B. G., Smith, A., Chen, C. M., and Silva, A. J. (1995). The α-Ca²⁺/calmodulin kinase II: a bidirectional modulator of presynaptic plasticity. Neuron 14, 591–597.

PubMed Abstract | Google Scholar

Chin, D., and Means, A. R. (2000). Calmodulin: a prototypical calcium sensor. Trends Cell Biol. 10, 322–328. doi: 10.1016/S0962-8924(00)01800-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Clapham, D. E. (2007). Calcium signaling. Cell 131, 1047–1058. doi: 10.1016/j.cell.2007.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, S. E., and Maler, L. (2017). Feedback synthesizes neural codes for motion. Curr. Biol. 27, 1356–1361. doi: 10.1016/j.cub.2017.03.068

PubMed Abstract | CrossRef Full Text | Google Scholar

Coultrap, S. J., and Bayer, K. U. (2012). CaMKII regulation in information processing and storage. Trends Neurosci. 35, 607–618. doi: 10.1016/j.tins.2012.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Csicsvari, J., Hirase, H., Mamiya, A., and Buzsáki, G. (2000). Ensemble patterns of hippocampal CA3-CA1 neurons during sharp wave-associated population events. Neuron. 28, 585–594. doi: 10.1016/S0896-6273(00)00135-5

PubMed Abstract | CrossRef Full Text | Google Scholar

de Jong, A. P. H., and Fioravante, D. (2014).Translating neuronal activity at the synapse: presynaptic calcium sensors in short-term plasticity. Front. Cell. Neurosci. 8:356. doi: 10.3389/fncel.2014.00356

CrossRef Full Text

De Koninck, P., and Schulman, H. (1998). Sensitivity of CaM kinase II to the frequency of Ca²⁺ oscillations. Science 279, 227–230.

PubMed Abstract | Google Scholar

Debanne, D., Bialowas, A., and Rama, S. (2013). What are the mechanisms for analogue and digital signalling in the brain? Nat. Rev. Neurosci. 14, 63–69. doi: 10.1038/nrn3361

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobrunz, L. E., and Stevens, C. F. (1997). Heterogeneity of release probability, facilitation, and depletion at central synapses. Neuron 18, 995–1008. doi: 10.1016/S0896-6273(00)80338-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupont, G., Houart, G., and De Koninck, P. (2003). Sensitivity of CaM kinase II to the frequency of Ca2+ oscillations: a simple model. Cell Calcium. 34, 485–497. doi: 10.1016/S0143-4160(03)00152-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Edelstein-Keshet, L. (2005). Mathematical Models in Biology. Philadelphia, PA: SIAM.

Google Scholar

Ferrell, J. E. (1997). How responses get more switch like as you move dawn a protein kinase cascade. Trends Biochem. Sci. 22, 288–289. doi: 10.1016/S0968-0004(97)82217-7

CrossRef Full Text | Google Scholar

Ferrell, J. E. Jr., and Ha, S. H. (2014). Ultrasensitivity part III: cascades, bistable switches, and oscillators. Trends Biochem. Sci. 39, 612–618. doi: 10.1016/j.tibs.2014.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Freedman, D., and Diaconis, P. (1981). On the histogram as a density estimator-L2 theory. Zeitschrift Wahrscheinlichkeitstheorie Verwandte Gebiete 57, 453–476. doi: 10.1007/BF01025868

CrossRef Full Text | Google Scholar

Graupner, M., and Brunel, N. (2007). STDP in a bistable synapse model based on CaMKII and associated signaling pathways. PLoS Comput. Biol. 3, 2299–2323. doi: 10.1371/journal.pcbi.0030221

PubMed Abstract | CrossRef Full Text | Google Scholar

Graupner, M., and Brunel, N. (2012). Calcium-based plasticity model explains sensitivity of synaptic changes to spike pattern, rate, and dendritic location. Proc. Natl. Acad. Sci. U.S.A. 109, 3991–3996. doi: 10.1073/pnas.110935910

PubMed Abstract | CrossRef Full Text | Google Scholar

Heltberg, M., Kellogg, R. A., Krishna, S., Tay, S., and Jensen, M. H. (2016). Noise induces hopping between NF-kappa B entrainment modes. Cell Systems 3, 532–539.e3. doi: 10.1016/j.cels.2016.11.014

CrossRef Full Text | Google Scholar

Hinds, H. L., Goussakov, I., Nakazawa, K., Tonegawa, S., and Bolshakov, V. Y. (2003). Essential function of alpha-calcium/calmodulin-dependent protein kinase II in neurotransmitter release at a glutamatergic central synapse. Proc. Natl. Acad. Sci. U.S.A. 100, 4275–4280. doi: 10.1073/pnas.0530202100

PubMed Abstract | CrossRef Full Text | Google Scholar

Hojjati, M. R., van Woerden, G. M., Tyler, W. J., Giese, K. P., Silva, A. J., Pozzo-Miller, L., et al. (2007). Kinase activity is not required for alpha CaMKII-dependent presynaptic plasticity at CA3-CA1 synapses. Nat. Neurosci. 10, 1125–1127. doi: 10.1038/nn1946

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, T. (1987). A thousand and one protein kinases. Cell 50, 823–829.

PubMed Abstract | Google Scholar

Jackman, S. L., Turecek, J., Belinsky, J. E., and Regehr, W. G. (2016). The calcium sensor synaptotagmin 7 is required for synaptic facilitation. Nature 529, 88–91. doi: 10.1038/nature16507

PubMed Abstract | CrossRef Full Text | Google Scholar

Jolivet, R., Lewis, T. J., and Gerstner, W. (2004). Generalized integrate-and-fire models of neuronal activity approximate spike trains of a detailed model to a high degree of accuracy. J. Neurophysiol. 92, 959–976. doi: 10.1152/jn.00190.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Jolivet, R., Schürmann, F., Berger, T. K., Naud, R., Gerstner, W., Roth, A., et al. (2008). The quantitative single-neuron modeling competition. Biol. Cybernetics 99, 417–426. doi: 10.1007/s00422-008-0261-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, Y. M., Chen, X. I., Lin, X. D., and Tan, N. (2017). Mean first passage time and stochastic resonance in a transcriptional regulatory system with non-gaussian noise. Fluct Noise Lett. 16:1750007. doi: 10.1142/S0219477517500079

CrossRef Full Text | Google Scholar

Kotaleski, J. H., and Blackwell, K. T. (2010).Modelling the molecular mechanisms of synaptic plasticity using systems biology approaches. Nat. Rev. Neurosci. 11, 239–251. doi: 10.1038/nrn2807

CrossRef Full Text

Larkman, A., Hannay, T., Stratford, K., and Jack, J. (1992). Presynaptic release probability influences the locus of long-term potentiation. Nature 360, 70–73.

PubMed Abstract | Google Scholar

Lewis, J. J., Slack, M. W., and Wolpert, L. (1977). Thresholds in development. J. Theor. Biol. 65, 579–590. doi: 10.1016/0022-5193(77)90216-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Stefan, M. I., and Le Novere, N. (2012). Calcium input frequency, duration and amplitude differentially modulate the relative activation of calcineurin and CaMKI. PLoS ONE 7:e43810. doi: 10.1371/journal.pone.0043810

CrossRef Full Text | Google Scholar

Lisman, J., Schulman, H., and Cline, H. (2002).The molecular basis of CaMKII function in synaptic and behavioral memory. Nat. Neurosci. Rev. 3, 175–190. doi: 10.1038/nrn753

CrossRef Full Text

Lisman, J., Yasuda, R., and Raghavachari, S. (2012).Mechanisms of CaMKII action in long-term potentiation. Nat. Rev. Neurosci. 13, 169–182. doi: 10.1038/nrn3192

CrossRef Full Text

Meyer, T., Hanson, P. I., Stryer, L., and Schulman, H. (1992). Calmodulin trapping by calcium-calmodulin dependent protein-kinase. Science 256, 1199–1202.

PubMed Abstract | Google Scholar

Michalski, P. J. (2013).The delicate bistability of CaMKII. Biophys. J. 105, 794–806. doi: 10.1016/j.bpj.2013.06.038

CrossRef Full Text

Michalski, P. J. (2014).First demonstration of bistability in CaMKII, a memory-related kinase. Biophys. J. 106, 1233–1235. doi: 10.1016/j.bpj.2014.01.037

CrossRef Full Text

Mizuseki, K., Royer, S., Diba, K., and Buzsáki, G. (2012). Activity dynamics and behavioral correlates of CA3 and CA1 hippocampal pyramidal neurons. Hippocampus 22, 1659–1680. doi: 10.1002/hipo.22002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ninan, I., and Arancio, O. (2004). Presynaptic CaMKII is necessary for synaptic plasticity in cultured hippocampal neurons. Neuron 42, 129–141. doi: 10.1016/S0896-6273(04)00143-6

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Leary, T., Williams, A. H., Caplan, J. S., and Marder, E. (2013). Correlations in ion channel expression emerge from homeostatic tuning rules. Proc. Natl. Acad. Sci. U.S.A. 110, E2645–E2654. doi: 10.1073/pnas.1309966110

PubMed Abstract | CrossRef Full Text | Google Scholar

Oswald, A. M., Lewis, J. E., and Maler, L. (2002). Dynamically interacting processes underlie synaptic plasticity in a feedback pathway. J. Neurophysiol. 87, 2450–2463. doi: 10.1152/jn.00711.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pang, Z. P., Cao, P., Xu, W., and Südhof, T. C. (2010). Calmodulin controls synaptic strength via presynaptic activation of calmodulin kinase II. J. Neurosci. 30, 4132–4142. doi: 10.1523/JNEUROSCI.3129-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, T. M., Schilstra, M. J., and Steuber, V. (2012). “The effective calcium/calmodulin concentration determines the sensitivity of CaMKII to the frequency of calcium oscillations,” in International Conference on Information Processing in Cells and Tissues (Berlin; Heidelberg: Springer). doi: 10.1007/978-3-642-28792-3_17

CrossRef Full Text | Google Scholar

Regehr, W. G. (2012). Short-term presynaptic plasticity. Cold Spring Harbor Perspect. Biol. 4:a005702. doi: 10.1101/cshperspect.a005702

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossetti, T., Banerjee, S., Kim, C., Leubner, M., Lamar, C., Gupta, P., et al. (2017). Memory erasure experiments indicate a critical role of CaMKII in memory storage. Neuron 96, 207–216 e2. doi: 10.1016/j.neuron.2017.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabatini, B. L., Oertner, T. G., and Svoboda, K. (2002). The life cycle of Ca2+ ions in dendritic spines. Neuron 33, 439–452. doi: 10.1016/S0896-6273(02)00573-1

CrossRef Full Text | Google Scholar

Sharma, G., and Vijayaraghavan, S. (2003). Modulation of presynaptic store calcium induces release of glutamate and postsynaptic firing. Neuron 38, 929–939. doi: 10.1016/S0896-6273(03)00322-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinha, S. R., Wu, L. G., and Saggau, P. (1997). Presynaptic calcium dynamics and transmitter release evoked by single action potentials at mammalian central synapses. Biophys. J. 72, 637–651.

PubMed Abstract | Google Scholar

Smedler, E., and Uhlen, P. (2014).Frequency decoding of calcium oscillations. Biochim. Biophys. Acta 1840, 964–969. doi: 10.1016/j.bbagen.2013.11.015

CrossRef Full Text

Smolen, P., Baxter, D. A., and Byrne, J. H. (1998). Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory systems. Am. J. Physiol. Cell Physiol. 274, C531–C542. doi: 10.1152/ajpcell.1998.274.2.C531

PubMed Abstract | CrossRef Full Text | Google Scholar

Stratton, M., Lee, I. H., Bhattacharyya, M., Christensen, S. M., Chao, L. H., Schulman, H., et al. (2014). Activation-triggered subunit exchange between CaMKII holoenzymes facilitates the spread of kinase activity. Elife 3:e01610. doi: 10.7554/eLife.01610

CrossRef Full Text | Google Scholar

Tank, D. W., Regehr, W. G., and Delaney, K. R. (1995). A quantitative analysis of presynaptic calcium dynamics that contribute to short term enhancement. J. Neurosci. 15, 7940–7952. doi: 10.1523/JNEUROSCI.15-12-07940.1995

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, L., Xie, Q., Ding, Y. H., Li, S. T., Peng, S., Zhang, Y. P., et al. (2013). CAMKII and calcineurin regulate the lifespan of Caenorhabditis elegans through the FOXO transcription factor DAF-16. Elife 2:e00518. doi: 10.7554/eLife.00518

PubMed Abstract | CrossRef Full Text | Google Scholar

Thanawala, M. S., and Regehr, W. G. (2013). Presynaptic calcium influx controls neurotransmitter release in part by regulating the effective size of the readily releasable pool. J. Neurosci. 33, 4625–4633. doi: 10.1523/JNEUROSCI.4031-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, G. M., and Huganir, R. L. (2004). MAPK cascade signalling and synaptic plasticity. Nat. Rev. Neurosci. 5, 173–183. doi: 10.1038/nrn1346

PubMed Abstract | CrossRef Full Text | Google Scholar

Urakubo, H., Sato, M., Ishii, S., and Kuroda, S. (2014). In vitro reconstitution of a CaMKII memory switch by an NMDA receptor-derived peptide. Biophys. J. 106, 1414–1420. doi: 10.1016/j.bpj.2014.01.026

PubMed Abstract | CrossRef Full Text | Google Scholar

von Wegner, F., Wieder, N., and Fink, R. H. (2014). Microdomain calcium fluctuations as a colored noise process. Front. Genet. 5:376. doi: 10.3389/fgene.2014.00376

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., and Maler, L. (1998). Differential role of Ca²⁺/calmodulin-dependent kinases in posttetanic potentiation at input selective glutamatergic pathways. Proc. Natl. Acad. Sci. U.S.A. 95, 7133–7138. doi: 10.1073/pnas.95.12.7133

CrossRef Full Text | Google Scholar

Wang, Z. W. (2008). Regulation of synaptic transmission by presynaptic CaMKII and BK channels. Mol. Neurobiol. 38, 153–166. doi: 10.1007/s12035-008-8039-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberg, S. H., and Smith, G. D. (2014). The influence of Ca2+ buffers on free [Ca2+] fluctuations and the effective volume of Ca2+ microdomains. Biophys. J. 106, 2693–2709. doi: 10.1016/j.bpj.2014.04.045

CrossRef Full Text | Google Scholar

Wen, Z. X., Guirland, C., and Ming, G. L. (2004). A CaMKII/calcineurin switch controls the direction of Ca2+-dependent growth cone guidance. Neuron 43, 835–846. doi: 10.1016/j.neuron.2004.08.037

CrossRef Full Text | Google Scholar

Xiong, W., and Ferrell, J. E. (2003). A positive-feedback-based bistable 'memory module' that governs a cell fate decision. Nature 426, 460–465. doi: 10.1038/nature02089

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhabotinsky, A. M. (2000). Bistability in the Ca2+/calmodulin-dependent protein kinase-phosphatase system. Biophys. J. 79, 2211–2221. doi: 10.1016/S0006-3495(00)76469-1

CrossRef Full Text | Google Scholar

Zheng, X. D., Yang, X. Q., and Tao, Y. (2011). Bistability, probability transition rate and first-passage time in an autoactivating positive-feedback loop. PLoS ONE 6:e17104. doi: 10.1371/journal.pone.0017104

PubMed Abstract | CrossRef Full Text | Google Scholar