A multi-state model of the CaMKII dodecamer suggests a role for calmodulin in maintenance of autophosphorylation

Ca2+/calmodulin-dependent protein kinase II (CaMKII) accounts for up to 2 percent of all brain protein and is essential to memory function. CaMKII activity is known to regulate dynamic shifts in the size and signaling strength of neuronal connections, a process known as synaptic plasticity. Increasingly, computational models are used to explore synaptic plasticity and the mechanisms regulating CaMKII activity. Conventional modeling approaches may exclude biophysical detail due to the impractical number of state combinations that arise when explicitly monitoring the conformational changes, ligand binding, and phosphorylation events that occur on each of the CaMKII holoenzyme’s subunits. To manage the combinatorial explosion without necessitating bias or loss in biological accuracy, we use a specialized syntax in the software MCell to create a rule-based model of a twelve-subunit CaMKII holoenzyme. Here we validate the rule-based model against previous experimental measures of CaMKII activity and investigate molecular mechanisms of CaMKII regulation. Specifically, we explore how Ca2+/CaM-binding may both stabilize CaMKII subunit activation and regulate maintenance of CaMKII autophosphorylation. Noting that Ca2+/CaM and protein phosphatases bind CaMKII at nearby or overlapping sites, we compare model scenarios in which Ca2+/CaM and protein phosphatase do or do not structurally exclude each other’s binding to CaMKII. Our results suggest a functional mechanism for the so-called “CaM trapping” phenomenon, wherein Ca2+/CaM may structurally exclude phosphatase binding and thereby prolong CaMKII autophosphorylation. We conclude that structural protection of autophosphorylated CaMKII by Ca2+/CaM may be an important mechanism for regulation of synaptic plasticity.


Introduction
CaMKII is a protein of interest because of its crucial role in synaptic plasticity [1][2][3][4][5]. In the hippocampus, synaptic plasticity in the post-synapse occurs within mushroom-shaped protrusions called dendritic spines [6]. Synaptic plasticity is dependent on calcium ion (Ca 2+ ) flux through N-methyl-D-aspartate receptors (NMDARs) located on the dendritic spines of the post-synaptic neuron [7]. Depending on the magnitude, frequency, and location of Ca 2+ flux, synaptic plasticity may produce increases or decreases (or neither) in synaptic strength [8,9]. Large, higher-frequency Ca 2+ spikes can induce an enduring up-regulation of synaptic strength, called long-term potentiation (LTP); while weak, lower-frequency Ca 2+ spikes can induce an enduring down-regulation of synaptic strength, called long-term depression (LTD) [9,10]. Whether Ca 2+ spikes induce LTP or LTD depends on relative activation of intracellular protein signaling networks. When Ca 2+ first enters the dendritic spine, it interacts with a variety of buffer and sensor proteins, chiefly calmodulin (CaM), which has many protein targets in the spine, including CaMKII [5,11,12].
The CaMKII holoenzyme contains at least twelve subunits [13][14][15][16] arranged as two rings of six. As shown in Fig 1, each CaMKII subunit features an N-terminal kinase domain and C-terminal hub domain [17]. Between the kinase and hub domains is a flexible regulatory domain which lends to the subunit a wide range of movement away from the holoenzyme's central hub. A crystal structure of human alpha-CaMKII expressed in E. coli published by Chao et al. (2011) shows CaMKII subunits as able to rapidly and stochastically pivot between a "docked" and "undocked" conformation, seemingly mediated by residues on the kinase domain's activation loop and a spur structure on the hub domain (see Fig 3C in [17]), such that a docked subunit may be inaccessible to CaM binding. In contrast, a more recent work using electron microscopy with rat alpha-CaMKII expressed in Sf9 cells suggests that less than 3 percent of subunits exhibit a compact (or docked) conformation [18]. Given the uncertainty in the field, we include subunit docking and undocking in our model, allowing for future exploration of this possible subunit functionality. In addition to docking and undocking, each subunit can be in an "inactive" conformation when the regulatory domain is bound to the kinase domain ( Fig  1B), or an "active" conformation when this binding is disrupted by thermodynamic effects. In this work, we describe CaMKII subunit activation as highly transient until stabilized by a fullybound Ca 2+ /CaM or phosphorylation at Thr-286 [17,19]. In the active conformation the catalytic domain of a subunit is able to bind and phosphorylate enzymatic substrates. A subunit may spontaneously return to an inactive conformation in the absence of Ca 2+ /CaM or phosphorylation at Thr-286 [19].
CaMKII can become enzymatically active in the absence of Ca 2+ /CaM-binding via phosphorylation at Thr-286, which is required for LTP [3,20]. Importantly, this phenomenon is an autophosphorylation: it is thought to occur when an active subunit phosphorylates neighboring subunits within the same holoenzyme [21,22]. Autophosphorylation at Thr-286 ("pThr-286") is thought to provide structural stability to a subunit's active conformation (reviewed in [23]) [24]. Because CaMKII plays a key role in the induction of LTP, and ultimately learning and memory (reviewed in [4,8]), we seek to better understand the biochemical regulation of CaMKII activation and autophosphorylation via computational modeling.
To characterize the spatiotemporal regulation of CaMKII, experimental studies are increasingly complemented by computational models [15,17,[25][26][27][28]. Computational models of The "inactive" CaMKII subunit (derived from a crystal structure found on the Protein Data Bank, or PDB, entry code 3SOA) in which the regulatory domain (yellow) is closely associated with the kinase domain (blue). (C) A schematic of the "active" CaMKII subunit. The regulatory domain (yellow) swings away from the kinase domain (blue). This schematic was generated by manually modifying PDB entry 3SOA to illustrate how the regulatory domain may be available for Ca 2 Ca 2+ -dependent signaling implicate competition, binding kinetics, feedback loops, and spatial effects in regulating enzyme activation [7,12,24,29,30]. However, fully characterizing these and other mechanisms of CaMKII regulation is impeded by the challenge of accurately portraying the CaMKII holoenzyme. As described by previous work, combinatorial explosion can occur when modeling CaMKII (and similar biomolecules) activation because the protein exhibits a large number of functionally significant and not necessarily inter-dependent states [24,26,[31][32][33]. The large number of possible states of CaMKII can neither be explicitly specified nor efficiently evaluated with conventional mass action-based methods. Indeed, for just one CaMKII hexamer ring, we estimate a state space of~32 billion states, and for the full dodecamer approximately 10 20 possible states (See Text A in S1 Appendix). The numbers of possible CaMKII states far exceeds the number of CaMKII molecules in a dendritic spine, suggesting that some states rarely occur and thus likely contribute little to protein function. Previous models leverage this observation to reduce the model state space and provide valuable insight to CaMKII binding and autophosphorylation dynamics [24,[33][34][35][36]. However, for CaMKII it remains unclear which states functionally participate in synaptic plasticity. Reduced models can inadvertently obscure key mechanisms regulating CaMKII activation and autophosphorylation. To elucidate complex regulatory mechanisms, it may be necessary for models to provide for all possible states ab initio.
In this work, we use rule-based model specification and particle-based rule evaluation methods to overcome combinatorial explosion [26,32,37]. Rules are conditions, based primarily on experimental observations, that prescribe when an implicitly-defined reaction may occur. At a given iteration, only states that matter for the execution of a particular rule are explicitly declared. States that do not matter to a particular rule can be omitted, a principle that has been paraphrased as "don't care, don't write" [38]. We use rule-and particle-based methods within the spatial-stochastic software MCell 3.3 [30,39] to present a comprehensive multistate model of the CaMKII dodecamer. Other simulation platforms can also overcome combinatorial explosion through rule-based model specification (e.g. BioNetGen [40]) or networkfree approaches (e.g. NFsim [41]). Unlike other platforms, MCell 3.3 provides both spatial-stochastic and rule-based modeling, although multi-state molecules in MCell 3.3 cannot diffuse. We use MCell 3.3 in anticipation of future MCell versions accounting for multi-state molecule diffusion, and to eventually build on simulations with physiological dendritic spine geometries such as those by   [42].
Here, we validate this rule-based MCell model of CaMKII regulation against current descriptions of the Ca 2+ frequency-dependence of CaMKII activation. By varying the rules and model parameter values we can simulate different experimental manipulations of CaMKII interaction with Ca 2+ /CaM and phosphatase and thereby explore various mechanisms regulating CaMKII activity. In particular, we show that Ca 2+ /CaM is important not only for regulating activation of CaMKII but may also contribute to the maintenance of CaMKII phosphorylation at Thr-286. We hypothesize that by limiting access of phosphatases to CaM-KII Thr-286 (perhaps by steric hindrance), Ca 2+ /CaM may prolong the lifetime of the autophosphorylated state.
CaM and PP are modeled in MCell as conventional cytosolic molecules. Initially, CaM is modeled as having one of two states: un-bound apo-CaM which does not bind CaMKII, and fully-saturated Ca 2+ /CaM (four Ca 2+ bound to CaM) which does bind CaMKII similar to previous studies [29,46]. Notably, we and others have described the importance of sub-saturated Ca 2+ /CaM states with fewer than 4 Ca 2+ [12,24,33,[47][48][49]. Thus, we also explore the dynamics of Ca 2+ -CaM binding and the binding of sub-saturated Ca 2+ /CaM to CaMKII. However, accounting for sub-saturated Ca 2+ /CaM would here require a multi-state representation, and because multi-state molecules cannot diffuse in MCell 3.3, we simplify our Ca 2+ /CaM model to allow CaM to diffuse and interact with a (non-diffusing) multi-state representation of CaM-KII. PP is modeled as single-state protein that is constitutively active and able to bind autophosphorylated CaMKII subunits. Our representation of constitutively active PP is consistent with previous models such as that by Lisman and Zhabotinsky (2001) [50].
CaMKII is modeled as a multi-subunit, multi-state complex, defined using a specialized model syntax for complex molecules (COMPLEX_MOLECULE) in MCell 3.3 [51]. This syntax allows for explicit representation of individual CaMKII dodecamers with distinguishable subunits. As shown in Fig 2, the holoenzyme is arranged as two directly-apposed, radially-symmetric rings each with six subunits. Each subunit can independently undergo state changes (e.g. subunit docking/undocking, CaM binding, phosphorylation at Thr-298, phosphorylation at Thr-306, binding to PP, see Fig 2). Each of these subunit state changes are represented as five "flags", each standing for a particular state that each CaMKII subunit can adopt. Note that all states are not mutually exclusive (i.e. a subunit can be phosphorylated at Thr-286 and bound to CaM simultaneously). Flags are used in rule evaluation, which occurs at each time step and for each individual subunit. That is, MCell repeatedly evaluates model rules against a given subunit's flags (and those of the neighboring subunits) to determine which state transitions a subunit undertakes at each time step. In the following sub-sections, we describe all CaMKII model flags, the state transitions that apply to each flag, the conditions and rate parameters for each state transition, and related model assumptions. In Fig 2, we visually convey how CaMKII subunits transition between states according to the model's rules. We summarize the state transition rules and rate parameter values in Table A in S1 Appendix.
Flag 1: Subunit docking. Docking is a binary flag that describes CaMKII kinase subunits as either "docked" or "undocked" to the CaMKII central hub. Subunits are instantiated in a docked state but may undergo numerous transitions between docked and undocked over the course of a simulation. At each time step, we assess a rule governing the subunit's transition from a docked to undocked state (see Table A in S1 Appendix). If this rule is satisfied, meaning that the subunit's docking flag is verified as "docked", then the transition to "undocked" is considered. Similarly, we assess a separate rule governing a transition from an undocked to docked state, which requires that the subunit not be bound to CaM and not phosphorylated at Thr-306 [17].
Subunit docking follows the structural model of Chao et al., who showed that a subunit cannot bind CaM as long as the subunit is in a compact conformation, docked to its central hub [17]. Docking implies a two-step process in which the subunit must first un-dock before subsequent CaM-binding, which accounts for the reported difference in binding rate for CaM to CaMKII-derived peptide (1 × 10 8 M -1 s -1 [52]) and for CaM to full-length CaMKII-T286A (1.8 × 10 6 M -1 s -1 [53]). Taking the ratio of these two rates gives an equilibrium constant for docking of 0.018, which is consistent with estimates by Chao et al., who assumed K docking to fall between 0.01 and 100 [17]. With this equilibrium constant, we estimate kinetic rates for docking and undocking. For this estimation, we first note that subunit docking involves a structural conformation change on a relatively large scale. Referring to a separate, and notably smaller-scale, conformational change in our model, in which CaM quickly transitions from an initially-to fully-bound state (see Flag 3: CaM Binding), we assume the docked-to-undocked transition to proceed at an order of magnitude slower. We therefore arrive at an assumed rate for k dock of 35 s -1 . In turn, this gives an undocking rate k undock = k dock × K docking of 0.63 s -1 , which lies within the range of 0.01 s -1 and 100 s -1 for k undock assumed by Chao et al.
Flag 2: Subunit activation. The activation flag describes CaMKII subunits as either "active" or "inactive". An inactive subunit has no catalytic activity because the regulatory domain is bound to the subunit's catalytic site. Conversely, an active subunit has catalytic activity because the regulatory domain's inhibition of the kinase domain is disrupted. When a subunit is active, Ca 2+ /CaM and/or other proteins may access and bind CaMKII. In our model, the transition reaction from inactive to active involves no explicit rules (but rather occurs unconditionally and as governed by rates described below). In contrast, two rules inform the conditions for subunit inactivation: that the subunit is 1) not fully-bound to CaM, and 2) not phosphorylated at Thr-286.
To assign rate parameters for this flag, we first note that subunits can fluctuate between inactive and active states rapidly in the absence of Ca 2+ /CaM (on the order of hundreds of nanoseconds) [19,54]. Noting this, we set the rate parameter for subunit inactivation at 1 × 10 7 s -1 . Further, Stefan et al. determined that the activation probability (in the absence of CaM and phosphorylation) is 0.002, leading us to set our activation rate parameter to 2 × 10 4 s -1 [31]. Thus, we arrive at a model in which CaMKII subunit activation is unstable until stabilized by CaM-binding or autophosphorylation.
Flag 3: CaM binding. CaM binding is a ternary flag meaning that each CaMKII subunit displays one of three states, where CaM may be "unbound", "initially-bound", or "fully bound". We emphasize that until we present a so-called "9-state-1-step" model of CaM-CaM-KII (see the end of Results), in this paper "CaM" generally refers to fully-saturated Ca 2+ /CaM (or CaM 4 ) bound to CaMKII. Our model adapts previous work by Stefan et al. (2012) to describe CaM-binding to CaMKII as a two-step process [31]. First, CaM binds to the regulatory domain of a CaMKII subunit (residues 298-312), resulting in a low-affinity "initially bound" CaMKII state, which is compatible with both the inactive (closed) and active (open) subunit conformation. Second, if the initially bound CaMKII opens it may transition to a "fully bound" state that describes the complete, higher-affinity interaction between CaM and CaMKII along residues 291-312 (see Fig 5 in [31]). The transition from an unbound to initially bound state requires the subunit must be: 1) undocked, 2) not PP-bound, and 3) un-phosphorylated at Thr-306. The transition reaction from initially bound to a fully bound state requires just that the subunit already be active/open. Dissociation of CaM from a fully bound CaM--CaMKII state proceeds through the initially bound state before becoming completely unbound from CaMKII.
In order to determine the parameters governing initial binding of CaM to CaMKII, we use data on CaM binding to CaMKII-derived peptides, rather than full-length CaMKII. This is done to separate the intrinsic binding constants from the parameters governing subunit activation/inactivation and docking/undocking. We use the microscopic k on for CaM binding to CaMKII measured by [52], using a CaMKII peptide and fluorescently labeled DA-CaM, as 1 × 10 8 M -1 s -1 . For the K D governing initial CaM binding, we use the K D reported by Tse et al. for CaM binding to a low-affinity peptide (CaMKII residues 300-312), which is 5.9 × 10 −6 M [55]. From these two parameters, we can compute the dissociation rate of initially-bound CaM from CaMKII: k off_CaM_ini = K d_CaM_ini × k on_CaM = 590 s -1 . While this rate may seem fast, we emphasize that in our model CaM dissociation happens in two steps, with the transition from fully-bound to initially-bound CaM which we discuss next.
In order to determine the parameters governing the transition from initially-bound to fully-bound CaM to CaMKII, we note that this transition involves a structural compaction of the CaM molecule, which has been measured using fluorescent labels [52,53]. Using fluorescent labels to analyze the structural compaction of CaM is convenient in its exclusion of effects due to conformational changes within CaMKII subunits or the CaMKII holoenzyme. Thus, we use these measurements as a proxy for CaM binding to a CaMKII peptide and to estimate parameters governing the transition between initially-and fully-bound CaM-CaMKII. Taken from experimental measurements by Torok et al., we identify a transition rate from initiallybound to fully-bound CaM-CaMKII (compaction of CaM) of 350 s -1 and from fully-bound back to initially-bound CaM-CaMKII (de-compaction of CaM) of 4 × 10 −3 s -1 [52]. This means that, in the absence of obstructions to binding, the likelihood of a bound CaM molecule being in the initial binding state (rather than the fully bound state) is 4 × 10 −3 / 350 = 1.1 × 10 −5 . This is consistent with a probability of CaM being bound to the high-affinity site of 0.99999 which was derived by Stefan et al. (2012) [31].
Flag 4: Phosphorylation at Thr-286. Phosphorylation at the residue Thr-286 is a ternary flag that describes this site as either "un-phosphorylated (uThr-286)", "phosphorylated (pThr-286)", or "phosphatase-bound". We specify four rules to govern the reaction that transitions a subunit from uThr-286 to pThr-286: the subunit 1) be uThr-286, 2) be undocked, 3) be active, and 4) have an active and undocked neighbor subunit in the same holoenzyme ring. Notably, we do not require that a subunit be CaM-bound for autophosphorylation to occur; however, because an un-bound CaMKII subunit's activation is highly transient, we find that CaM-binding is effectively required for autophosphorylation at Thr-286 (shown in Model Validation). Similarly, because pThr-306 (discussed in Flag 5) prohibits CaM-binding, we find that autophosphorylation at Thr-286 effectively requires uThr-306, though we do not explicitly state this rule. The neighboring subunit's activation flag is considered because autophosphorylation is facilitated by its catalytic site. Our model only considers the counter-clockwise neighbor subunit because, in the absence of experimental observations to the contrary, we assume that steric effects cause autophosphorylation to occur in only one direction about a CaM-KII ring, similar to previous work [56,57]. The rate of autophosphorylation, 1 s -1 , at Thr-286 is taken from an earlier study of CaMKII autophosphorylation in the presence of CaM [47].
De-phosphorylation of pThr-286 is facilitated by binding and enzymatic activity of protein phosphatases PP1 and PP2A, here referred to generally as PP [43,44]. Two rules govern PP binding to a CaMKII subunit (the transition from pThr-286 to a phosphatase-bound state): that the subunit be 1) phosphorylated at Thr-286 (pThr-286) and 2) not be bound to CaM. It has been shown that a majority of autophosphorylated CaMKII in the PSD is dephosphorylated by PP1 [58,59]; while in brain extracts autophosphorylated CaMKII is mostly dephosphorylated by PP2A [43]. The requirement that CaM be unbound from CaMKII in order for PP to bind to CaMKII is motivated by the observation that simultaneous binding of CaM and PP to the CaMKII regulatory domain may be mutually exclusive due to steric hindrance. CaM, having molecular weight 18 kDa, binds to the CaMKII regulatory domain around residues 290-309 [56,60,61], which is at least 4 residues, and at most 23 residues away from Thr-286 (again, see also Fig 5 in [31]). To the best of our knowledge, the peptide binding footprint of neither PP (PP1 nor PP2A) onto CaMKII is fully described. However, both PP1 and PP2A are widely known to target pThr-286 [58,59,62] and de-phosphorylate threonine residues nearby alpha helices in other substrates [63,64]. Additionally, the catalytic subunit of PP1 has a molecular weight of 37 kDa, which is nearly twice that of CaM and more than half that of a CaMKII subunit. Taken together, we hypothesize that the PP binding footprint likely overlaps with the CaM binding site, such that the presence of bound PP likely structurally excludes or impedes upon a subsequent binding of CaM to CaMKII. Similarly, the presence of bound Ca 2+ /CaM structurally would exclude coincident binding of PP. Starting with Text C in S1 Appendix, we further discuss the quantitative basis of this structural exclusion hypothesis in light of the crystal structure of the PP1-spinophilin interaction (PDB: 3EGG) [65]. In short, PP1 tends to bind substrates at a site >20Å from the PP1 active site. Thus, if the PP1 binding footprint does not actually contain T286, then the furthest likely CaMKII residue of PP1 binding (at least on the hub domain side of T286) is G301, well within the CaM binding footprint (Fig J in S1 Appendix). We examine the regulatory implications of this hypothesis by relaxing the rules of PP binding and requiring only that the subunit be pThr-286. The kinetic rates of PP binding CaMKII are based on values used by Zhabotinsky (2000), using the same catalytic rate of 2 s -1 and choosing values for the association (3x10 6 M -1 s -1 ) and dissociation (0.5 s -1 ) rate constants such that the resulting Michaelis constant (8.3 x 10 −7 M) falls in the range used by Zhabotinsky [46]. Flag 5: Phosphorylation at Thr-306. Phosphorylation at the residue Thr-306 is a binary flag that describes this site as either un-phosphorylated ("uThr-306") or phosphorylated ("pThr-306"). We model the transition from uThr-306 to pThr-306 using three rules: that that the subunit be 1) uThr-306, 2) active, and 3) un-bound by CaM. Our model uses a forward rate parameter 50-fold slower than that of phosphorylation at Thr-286, based on past experimental measurements [35,66]. Over the course of our simulation times, we observe very few pThr-306 transitions and therefore exclude the reverse transition reaction describing de-phosphorylation of pThr-306 into uThr-306. Notably, we do not require that a subunit be pThr-286 as a condition for becoming pThr-306. However, because a subunit's activation is highly transient when un-bound by CaM in our model, the probability of pThr-306 acquisition is quite small if the subunit is not already pThr-286.

Model validation
Stimulation frequency correlates with subunit activity. To validate our model, we assessed a variety of model outputs under various regimes of Ca 2+ /CaM stimulation. As a first assessment, we simulated a persistent Ca 2+ /CaM bolus (note that here we only use the fully saturated CaM 4 as the active form of Ca 2+ /CaM), similar to experiments by Bradshaw et al. (2002), who monitored CaMKII autophosphorylation over time [57]. In Fig 3 we simultaneously monitored the time-course concentration of CaMKII subunit flags indicating: initially-bound Ca 2+ /CaM, fully-bound Ca 2+ /CaM, active CaMKII, and pThr-286 (Fig 3). In the persistent, continuous presence of Ca 2+ /CaM, the concentration of subunits with initiallybound Ca 2+ /CaM (yellow trace) is noisy and consistently low, implying that initially-bound Ca 2+ /CaM seems rapidly to either dissociate or proceed to a fully-bound conformation. Fullybound Ca 2+ /CaM (red trace) subunit concentrations closely follow those of active CaMKII subunits (dark blue trace) over time, providing evidence that Ca 2+ /CaM stabilizes CaMKII activation. Indeed, because the difference in concentrations of fully-bound Ca 2+ /CaM and active CaMKII is always small, we observe that although unbound CaMKII may spontaneously activate, these activated subunits rapidly return to an inactive state and are not likely to progress to a phosphorylated (pThr-286) state. We next observe that the increase of CaMKII autophosphorylation at Thr-286 (cyan trace) over time is strongly associated with increases in the number of subunits that are fully-bound to Ca 2+ /CaM and active subunits (dark blue and red traces, respectively). This is consistent with previous work showing that Ca 2+ /CaM must be bound to CaMKII for pThr-286 to occur [56] and CaMKII Ca 2+ -independent activity is strongly associated CaMKII autophosphorylation at Thr-286 [17,53,67,68]. Furthermore, we observe in Fig 3A that [57]).
Next, we assessed model behavior under low-and high-frequency stimulating conditions. CaMKII activation and autophosphorylation at Thr-286 in response to 5Hz and 50Hz Ca 2+ /CaM is plotted in Fig 3B and 3C, respectively; 50 seeds were run for each condition, with 6 representative traces (transparent lines) and the average response (bold) plotted. As expected, the data showed significantly greater levels of CaMKII activation and autophosphorylation at 50Hz relative to 5Hz stimulation [12,20]. Indeed, our results in Fig 3C are comparable to results from Shifman et al. (2006), who observed much lower autophosphorylation at low Ca 2+ /CaM concentrations (less than 2 μM) than at high concentrations (see Fig 4D in [47]).
To further determine how stimulation frequency affects CaMKII activity, the model was stimulated continuously at frequencies ranging from 1Hz to 50 Hz. The concentrations of the various CaMKII states t = 20 seconds of simulation time are plotted as a function of frequency in Fig 3D. We observe a nearly linear correlation between both subunit activation (R 2 = 0.99) and pThr-286 (R 2 = 0.96). This is consistent with computational results from Chao et al., who developed a stochastic model that also yielded a linear relationship between pThr-286 and stimulation frequency for frequencies greater than 1 Hz [15]. Additionally, our results in Fig  3D show that the model elicits~5% (~1μM out of 18.24μM total) CaMKII pThr-286 in response to 10Hz stimulation (Fig 3D), which is in general agreement with experimental results given our relatively short pulse width of 10msec (see Fig 4A in [69]). In Fig A in S1 Appendix, we further simulate our model for 20sec at 4Hz Ca 2+ /CaM using a pulse width of 200msec, yielding~50% pThr-286, consistent with [69] (see Fig 3B in that paper). Taken together, these results (Fig 3 and Fig A in S1 Appendix) show that our model behaves as expected and is able to produce CaMKII activity and autophosphorylation behaviors similar to previous computational and experimental results.

Exploring CaMKII holoenzyme phosphorylation
A Thresholded response of CaMKII to Ca 2+ CaM. CaMKII has long been theorized to exhibit switch-like or bistable behavior, which could underlie the importance of pThr-286 to learning and memory formation [4,46,50,70,71]. However, experimental efforts have struggled to identify a bistability making CaMKII activity robust to the presence of phosphatases. Recently, Urakubo et al. used the chelator EGTA to control single pulses of Ca 2+ in a mixture of CaM, CaMKII, PP1, and NMDAR peptides, leading to what seemed to be the first direct observation of CaMKII bistability in the presence of NMDAR peptides only [72]. Noting this, we explored whether a spatial stochastic model of the CaMKII dodecamer would exhibit switch-like behavior for concentration parameters of Ca 2+ , CaM, CaMKII, and PP known to exist in hippocampal spines. We stimulated the model with a set of short Ca 2+ /CaM input pulses (which could also be reproducible in vitro). Importantly, we did not aim to identify true bistability because exploring the many combinations of Ca 2+ , kinase, and phosphatase concentrations was outside the scope of this paper. Instead we wondered if, by stimulating with brief pulses of Ca 2+ /CaM of variable duration, our model would exhibit switch-reminiscent pThr-286 behavior. Specifically, we hypothesized that there would exist a threshold of time of Ca 2+ /CaM below which pThr-286 was unachievable and above which pThr-286 was maintained. At stimulation durations above this threshold enough CaMKII autophosphorylation (pThr-286) would occur and self-propagate in balance with de-phosphorylation by PP.
In  (Fig 4A  and 4B, respectively). Interestingly, subunits stimulated by even the shortest pulses of 0.05 or 0.1 sec, appeared to sustain their activation for the simulation period (120 sec). However, these short-pulse (0.05-0.1 sec) stimulations rarely resulted in autophosphorylation (pThr-286, Fig  4B). Longer (0.2-0.5 sec) Ca 2+ /CaM pulses resulted in greater levels of subunit activation that started declining immediately after the Ca 2+ /CaM pulse ended (Fig 4A), but elicited pThr-286 levels that were generally sustained for the duration of a simulation (Fig 4B). Taken together, we found that CaMKII may be thresholded at a level of Ca 2+ /CaM exposure below which pThr-286 is unobserved and above which pThr-286 is achieved and subsequently sustained across several minutes even in the presence of phosphatase.
We also explored how this Ca 2+ /CaM threshold may depend on the number of directions by which subunits can autophosphorylate their neighbors. Note that in the results up to this point, autophosphorylation was limited to occurring in a single direction, or degree of freedom. That is, subunits could only autophosphorylate their adjacent neighbors in a counterclockwise fashion [17,56,57]. We therefore created alternative versions of our model in which autophosphorylation could occur with multiple degrees of freedom: two degrees of freedom in which intra-ring phosphorylation occurs in counter-and clock-wise directions, and three degrees of freedom in which both intra-ring and trans-ring phosphorylation occurs. We used these higher-degree of freedom models to monitor the rates of pThr-286 formation both in bulk and on an individual subunit basis. As expected, pThr-286 formation and intra-holoenzyme propagation rates increased with increasing degrees of freedom (see Figs E-G in S1 Appendix), though the differences would likely not be distinguishable by bench-top experimentation. In addition, the length of time in which consecutive neighboring subunits remained autophosphorylated also increased with increasing degrees of freedom (Fig G in S1  Appendix). This implied that subunits may be more frequently autophosphorylated on timeaverage with increasing degrees of freedom. We present this preliminary exploration of the implication of the degrees of freedom, or directionality of CaMKII holoenzyme autophosphorylation, in order to demonstrate various capabilities of the model. We note that current experimental techniques are not sensitive enough to be discriminated between the simulation results. More work needs to be done to better characterize the conditions or possibility of bidirectional autophosphorylation both experimentally and computationally. Future experimental and computational studies could perhaps explore the dependence of autophosphorylation on higher degrees of freedom.
CaM-dependent exclusion of PP1 binding stabilizes autophosphorylation. Fig 4 suggested a threshold of Ca 2+ /CaM activation beyond which CaMKII remains autophosphorylated, implying a balance between kinase and phosphatase activity. We wondered how a putative balance between CaMKII autophosphorylation and phosphatase activity might be regulated. In the previous experimental work by Urakubo et al., maximally-phosphorylated CaM-KII was maintained in the presence of PP1 and NMDAR peptide for as long as 8 hours (at 4˚C). In that work, addition of the kinase inhibitor K252a to phosphorylated CaMKII resulted in a steady decline in pThr-286 towards basal levels, suggesting that maintenance of pThr-286 over time was not due to low phosphatase activity, but rather a recovery of de-phosphorylated subunits back to a phosphorylated state. To recreate inhibition of kinase activity in our model, at time t = 30 sec we introduced a high concentration (18.2 μM) of K252a, enough to bind all CaMKII subunits in the model. K252a binding results in a blocked CaMKII state that cannot be autophosphorylated (see Flag 2 in Table A in S1 Appendix). Importantly, the blocked CaM-KII subunit can still be de-phosphorylated at pThr-286. In separate simulations we explored the effects of a phosphatase inhibitor, which was also introduced at t = 30sec. To simulate the introduction of a phosphatase inhibitor, we defined the catalytic rate of de-phosphorylation by PP1 (k cat PP1 ) as a time-dependent variable that assumed a value of zero at t = 30sec. This implementation of kinase and phosphatase inhibition preserved normal CaM and PP1 binding dynamics.
In these simulations the model was stimulated with 22.8μM CaM 4 (Ca 2+ /CaM) for 2sec and either no inhibition (control), kinase activity inhibition, or phosphatase inhibition was implemented at 30sec as described above (Fig 5). As expected, inhibiting phosphatase activity (green trace) caused kinase activity to dominate, resulting in a steady increase in pThr-286 compared to the control (black trace). Surprisingly, the kinase-inhibition (blue trace) showed little difference in pThr-286 concentration over time compared to the control. Instead of causing a greater reduction in pThr-286 over time due to phosphatase activity as was expected, inhibition of kinase activity resulted in almost no difference in pThr-286 levels even in the presence of phosphatase activity. We hypothesized that some other, non-enzymatic mechanism in our model was contributing to the maintenance of pThr-286.
In every simulation presented thus far, we assumed that CaM binding to the CaMKII regulatory domain sterically hinders PP binding to the regulatory domain, and vice-versa. This was implemented in the model via a rule that requires a subunit be unbound by CaM in order for PP to bind. To test the role of these exclusions, we created a second version of our model in which PP binding would be allowed regardless of the presence of bound CaM 4 , and CaM 4 binding would be allowed regardless of the presence of PP. In contrast to our original "exclusive" model, the "non-exclusive" model required only that a subunit be pThr-286 in order for PP binding to be allowable. In other words, the non-exclusive model allowed Ca 2+ /CaM and PP to bind CaMKII agnostically of each other. Aside from this rule adjustment, our exclusive and non-exclusive models utilized identical parameters (see Tables A and B in S1 Appendix). As in Fig 5, we selected a Ca 2+ /CaM bolus time of 2 sec. Again, we monitored both CaMKII activation (Fig 6A and 6B) and pThr-286 (Fig 6C and 6D) over 120 seconds of simulated time. Critically, both the exclusive and non-exclusive models were examined with high (purple trace) and low (orange trace) association rate parameter values for PP binding to CaM-KII. Increasing and decreasing the association rate of PP (k on PP is normally set to 3 μM -1 sec -1 ) to CaMKII by one order of magnitude accounted for parameter uncertainty and provided a magnified view of the signaling effects of CaM-mediated exclusion of PP binding. Our results suggested that CaM-dependent exclusion of PP is an important regulatory mechanism for maintaining CaMKII autophosphorylation levels. While the PP exclusion rule had little to no effect on the decay of CaMKII subunit de-activation ( Fig 6A and Fig 6B, both decay constants -0.0004), pThr-286 (Fig 6C and Fig 6D) was highly influenced by the PP exclusion rule. In the exclusive model (Fig 6C), pThr-286 levels were steady and stable (decay constant -0.0004) despite varying the PP association rate parameter by two orders of magnitude. In contrast, the non-exclusive model (Fig 6D) showed that for a high PP association rate, significant pThr-286 levels remained below 0.5 μM (Fig 6D low k on trace has decay constant -0.006). Moreover, for a low PP association rate, the non-exclusive model attained lower pThr-286 levels compared to the exclusive model, and the pThr-286 levels then declined at a faster rate. We further emphasize the magnitudes of pThr-286 were noticeably influenced by the PP exclusion rule. Upon removing PP exclusion, peak magnitudes of pThr-286 were reduced from 3.2μM (Fig 6C) to~2.3μM (Fig 6D). It seemed that in order to significantly activate and then maintain pThr-286 over longer time periods, CaMKII required a mechanism regulating phosphatase access, and a regulator of phosphatase access could be CaM itself.  ) with CaMKII is either increased (purple To reinforce our assertion that CaM-dependent structural exclusion of PP binding stabilizes pThr-286, we repeated simulations shown in Fig 5, but with our non-exclusive model. In Fig 6E, we stimulated our non-exclusive model with a 2sec pulse of Ca 2+ /CaM and then monitored pThr-286 over time. For these simulations, k on PP1 was restored to its standard value of 3 μM -1 sec -1 . As in Fig 5, in separate simulations we inhibited at t = 30sec either phosphatase activity, kinase activity, or neither (control). The control (grey trace) was reminiscent of results in Fig 6D, in which pThr-286 was achieved but then slowly declined on a steady yet noisy basis. Notably, all non-exclusive model variants were much noisier than their exclusive model counterparts in Fig 6E. Inhibiting phosphatase activity (light green trace) in the non-exclusive model again caused kinase activity to dominate and pThr-286 levels to generally increase over time, similarly to the exclusive model. In contrast to the exclusive model, inhibiting kinase activity (light blue trace) in the non-exclusive model rapidly and totally abolished pThr-286. It seemed that for the non-exclusive model, in which CaM and PP could bind simultaneously, inhibiting kinase activity caused phosphatase activity to dominate. Taken together, these results suggested that in addition to supporting CaMKII subunit activation, CaM also has a role in maintaining CaMKII activity by blocking phosphatase access and thereby slowing down dephosphorylation. CaM-dependent exclusion of PP1 binding may depend significantly on how we model Ca 2+ /CaM. Until this point, we have modeled Ca 2+ /CaM as "2-state-2-step", existing as either Ca 2+ -unbound apo-CaM or CaM 4 (2 states), which binds CaMKII in an initially-then fullybound conformation (2 steps). However, previous experimental and computational studies have determined that sub-saturated Ca 2+ /CaM, with fewer than four Ca 2+ bound, may significantly bind CaM-binding partners such as CaMKII [24,47]. Indeed, Pepke et al. [24] and others use a "9-state-1-step" model of Ca 2+ /CaM, which explicitly accounts for each mode of Ca 2+ -binding at the CaM N-and C-termini. Importantly, each of the nine Ca 2+ /CaM states in the Pepke model has unique binding kinetics for CaMKII. We emphasize that these 9-state binding kinetics, which were measured using wild-type CaMKII in vitro, are incompatible with our 2-step CaM-binding model. In other words, a 9-state-2-step CaM-CaMKII model is difficult to parameterize because the available 9-state parameter values inherently account for 2-step CaM-binding. Moreover, the 9-state-2-step model would likely require a multi-state, rule-based model of CaM. And problematically, MCell 3.3 prohibits diffusion for rule-based species. Still, it is important to consider whether sub-saturated Ca 2+ /CaM states might still be able to structurally exclude, or out-compete, PP1-CaMKII binding.
Although a 9-state-2-step model of CaM-CaMKII binding is currently impractical, a 9-state-1-step model of CaM-CaMKII binding is practical, at least to explore how sub-saturated Ca 2+ /CaM could exclude PP1-CaMKII binding. For the 9-state-1-step model, we again use MCell 3.3 to describe the multi-state CaMKII holoenzyme, but we modify three of the flags described earlier in this paper. First, we remove the subunit docking and activation flags to reduce model noise and ensure the 9-state CaM-binding parameters remain appropriate. Second, we modify the CaM-binding flag to allow all nine Ca 2+ /CaM states (including apo-CaM) to bind CaMKII subunits. Thus, whereas in the 2-state-2-step CaMKII model subunit activation is defined by subunit opening, in the 9-state-1-step model activation is defined by CaMbinding. The parameters and reaction network for Ca 2+ -CaM binding and CaM-CaMKII binding may be found in Pepke et al. (see Fig 2C in [24]), and also refer to Table C  Appendix. With our 9-state-1-step model (MCell code provided in the Purdue University Research Repository with DOI: 10.4231/MV0Z-8Z57, and parameters provided in S1 Appendix), we simulate using identical conditions to those used for Fig 6A-6D, namely with 2000 Ca 2+ (a CaM-saturating quantity, fully chelated at t = 2sec), 450 CaM particles, and 30 CaMKII holoenzymes. Our 9-state-1-step model results are shown in Fig 7, where we again show that CaMKII subunit activation and pThr-286 levels are maintained on significantly longer timescales when CaM and PP cannot bind CaMKII simultaneously (Fig 7A and 7D). Compared to Fig 6, activation and pThr-286 levels are both higher and output noise is reduced, likely due to the absence of the activation flag and the fact that in a 9-state model, sub-saturated Ca 2+ /CaM significantly bind CaMKII and could contribute to its autophosphorylation. These may also explain the lack of difference in the peak magnitudes of pThr-286 in Fig 7C and 7F (compare to the peak magnitudes of Fig 6B and 6D). Specifically, for the 9-state-1-step model in Fig 7  maximal phosphorylation is achievable within 2 sec of Ca 2+ stimulation, whereas 2 sec is not sufficient for maximal phosphorylation in the 2-state-2-step simulations (Fig 6). Finally, we observe that when only fully-saturated CaM 4 but no other Ca 2+ /CaM state is allowed to exclude PP (Fig 7B and 7E), model output is virtually identical to the fully non-exclusive case (Fig 7C and 7F). This suggests that sub-saturated Ca 2+ /CaM states may significantly contribute to PP1 exclusion. Indeed, this may be unsurprising given the affinity of states such as CaM 2C for CaMKII (7.4μM), which is only one order of magnitude larger than that of PP1 for CaMKII (0.166μM). Critically, though, in our 9-state-1-step model the CaM 2C -CaMKII affinity increases by 1000-fold when a CaMKII is subunit is pThr-286, because we explicitly account for CaM-trapping [73]. Therefore, sub-saturated Ca 2+ /CaM states are very likely to out-compete PP1 and prevent its binding to CaMKII, at least from a kinetics perspective.

Discussion
In this work, we use rule-and particle-based methods with the software MCell to model the complete CaMKII holoenzyme. Rule-based modeling allows us to account for and monitor multiple CaMKII states simultaneously without eliciting combinatorial explosion. By explicitly accounting for multiple CaMKII states, we can use this model to explore regulatory mechanisms such as the CaM-dependent maintenance of pThr-286 by structural exclusion of phosphatase binding to CaMKII.
Previous multi-state models of CaMKII exist but are different in focus and in scope from the present model. For example, our model is based on an earlier multi-state model by Stefan et al. (2012) [31] implemented in the particle-based stochastic simulator StochSim [74]. Stoch-Sim accounts for subunit topology (i.e. the user can specify whether a subunit is adjacent to another, and reactions can be neighbor-sensitive), but StochSim does not explicitly account for spatial information. MCell, as a spatial simulator, offers more possibilities to precisely account for spatial effects and to situate models in spatially realistic representations of cellular compartments. In addition, the model by Stefan et al. provides only for interactions between adjacent CaMKII molecules on the same hexamer ring and therefore models CaMKII as a hexamer, not a dodecamer. Similarly, another previous model of CaMKII by Michalski and Loew (2012) uses the softwares BioNetGen and VCell to offer an infinite subunit holoenzyme approximation (ISHA) of the CaMKII hexamer [75][76][77]. The ISHA asserts that under certain enzymatic assumptions, the output of a multi-state CaMKII model is independent of holoenzyme size when the number of subunits exceeds six. However, Michalski's ISHA model is most suitable for systems containing only one holoenzyme structure-dependent reaction such as the autophosphorylation at Thr-286. Additional reactions to describe actin binding [78] or subunit exchange [14,15] may invalidate Michalski's ISHA, whereas our model can in the future readily accommodate additional, holoenzyme structure-dependent phenomena. Finally, a more recent rule-based model of the CaMKII holoenzyme by Li and Holmes [26] offers a detailed representation of how CaM binds to Ca 2+ and subsequently activates CaMKII subunits, based on earlier results of CaM regulation [79]. Li and Holmes offer valuable and detailed insight into how CaM binding to CaMKII depends on Ca 2+ dynamics. While our model is less detailed in representing the regulation of CaM itself, our model is much more detailed in representing other aspects of CaMKII regulation, including multiple modes of CaM binding, conformational change, detailed holoenzyme structure, multiple phosphorylation sites, and dephosphorylation. We can in the future expand our MCell model to account for multiple holoenzyme structure-dependent phenomena and simultaneously incorporate this model into the broader Ca 2+ -dependent signaling network.
This work in-part demonstrates the value of MCell as a rule-based modeling framework. Rule-based modeling accommodates much larger state spaces than is possible using conventional systems of differential equations. Admittedly, not all models (including models of CaM-KII) require extensive state spaces, but rule-based modeling results can help justify the assumptions typically used to reduce a state space. For example, our model conditions yield, as shown in Fig 3A, negligible levels of initially-bound CaM compared to other states such as fully-bound CaM or pThr-286. Therefore, it might sometimes be appropriate to exclude an initially-bound CaM state from future implementations in frameworks for which combinatorial explosion is a concern. Aside from addressing combinatorial explosion, rule-based models are especially well-suited to discern otherwise concealed mechanisms, as exemplified by Di Camillo et al. who used rule-based models to identify a robustness-lending negative feedback mechanism in the insulin signaling pathway [49]. Furthermore, MCell describes CaMKII holoenzymes as discrete particles in space, which will lend realism to future spatial-stochastic models of Ca 2+ -dependent signaling networks in the dendritic spine, a compartment in which the Law of Mass Action is invalid [24]. This particle-based framework also allows for individual subunit monitoring, which works in conjunction with the Blender software plugin, CellBlender (see Movie A in S1 Appendix).
One of the results of this work is the identification of distinct levels of CaMKII activation and pThr-286 in response to distinct pulses of Ca 2+ /CaM stimulation. Distinct levels of CaM-KII activation could tune the selectivity of CaMKII for certain downstream binding targets such as AMPA receptors or the structural protein PSD-95. If stimulation-dependent tuning of CaMKII activation were observed, it would be reminiscent of other studies that have implicated feedback loops [29] and binding dynamics [24] as regulators of Ca 2+ -dependent enzyme activation. For example, a recent study suggests that competition is an emergent property that tunes the Ca 2+ frequency dependence of CaM binding to downstream targets, leading Ca 2 + /CaM to set distinct levels of calcineurin-and CaMKII-binding [12]. Similarly, CaMKII itself could preferentially select downstream binding partners as a function of its level of activation by Ca 2+ /CaM, possibly providing a mechanism by which CaMKII facilitates certain LTPrelated molecular events. Additionally, our observation of distinct levels of CaMKII activation and thresholded pThr-286 could be an indication of long-hypothesized switch-like behavior in synaptic plasticity [4,70]. If switch-like behavior in fact occurs, then pThr-286 is likely maintained by a balance in kinase and phosphatase activity.
While investigating a putative interplay in CaMKII kinase and PP phosphatase activity in maintaining pThr-286 levels, we may have identified a CaM-dependent mechanism that blocks PP binding to CaMKII. In a model that excludes simultaneous binding of CaM and PP to CaMKII, pThr-286 significantly increases upon phosphatase inhibition, yet in the same model kinase inhibition causes little change in pThr-286 over time (Fig 5). In contrast, a nonexclusive model that allows simultaneous binding of CaM and PP shows that introduction of a kinase inhibitor rapidly abolishes pThr-286. These results suggest that CaM-dependent exclusion of PP may provide a stabilizing mechanism. Additionally, we use our MCell-based implementation of the model to monitor transitions between multiple states of distinct subunits within holoenzymes (Fig 8 and Movie A in S1 Appendix).
The major advance of this paper is to present a model of the regulation of 12-subunit CaM-KII holoenzyme and its regulation by Ca 2+ /CaM and protein phosphatase. We assert that CaM-dependent exclusion of PP could provide a functional role for so-called "CaM trapping" [56] and possibly contribute to CaMKII bistability. Indeed, a model by Zhabotinsky (2000) explored CaMKII bistability, indicating that two stable states of pThr-286 would in-part require very high CaMKII concentrations, seemingly to bolster kinase activity in the system [46]. However, the Zhabotinsky model assumes that CaM and PP1 could bind CaMKII simultaneously, possibly exaggerating the ability of PP1 to de-phosphorylate at Thr-286. If PP1 binding were to be encumbered in the Zhabotinsky model, perhaps through CaM-dependent exclusion, then bistability might be achievable at lower CaMKII concentrations.
Previous studies have sought to explore the dependence of CaMKII de-phosphorylation on the presence of Ca 2+ /CaM. An experiment by Bradshaw et al. (2003) quantifies PP1-mediated de-phosphorylation rates of pThr-286 in vitro, in the presence or absence of free Ca 2+ (see Fig  4B in [80]). In both the presence and absence of free Ca 2+ at 0˚C, Bradshaw et al. observe that pThr-286 levels decrease over 30 min at similar rates, in conditions with limiting ATP such that CaMKII cannot re-phosphorylate. These results could be interpreted to suggest that PP1 activity is independent of Ca 2+ and in turn Ca 2+ /CaM-binding to CaMKII. Note that in the Bradshaw results, regardless of the presence of free Ca 2+ , CaMKII activation persists for at least tens of minutes. The persistence of CaMKII activity on relatively long timescales, even in the presence of phosphatase, is consistent with separate experimental results [22,57,59,70,73]. Our results suggest an alternative interpretation of these and Bradshaw's results. In our model, CaMKII activity persists for many minutes only when CaM excludes PP1. Indeed, the effect of CaM exclusion on the timescale of CaMKII activation is even more pronounced in our 9-state-1-step CaM-binding model results.
Our 9-state-1-step model suggests that following termination of Ca 2+ flux, sub-saturated Ca 2+ /CaM states may significantly contribute to PP1-exclusion (Fig K in S1 Appendix). First, this agreement between our 9-state-1-step and primary 2-state-2-step models is a testament to the model's robustness to uncertainty with parameterizing the docking and activation flags. Further, our 9-state-1-step results are also consistent with Bradshaw's results, as seen in Fig 9 of [57] (their supplementary material). Specifically, Bradshaw et al. explore in their supplement how CaMKII autophosphorylation levels equilibrate in the presence of various Ca 2+ concentrations. Importantly, many of the Ca 2+ concentrations are at levels that would not saturate CaM. Nonetheless, with sub-saturated Ca 2+ /CaM significant levels of CaMKII autophosphorylation are maintained even by 7.5 hours of incubation, even in the presence of PP1, indicating that sub-saturated Ca 2+ /CaM states could be playing a key role in maintenance of CaMKII autophosphorylation. In S1 Appendix, we show that the Ca 2+ /CaM state predominantly responsible for PP1-exclusion following termination of Ca 2+ stimulation in the 9-state-1-step model is apo-CaM. Apo-CaM remaining bound to CaMKII is consistent with results by Brown et al., who determine that when free Ca 2+ levels decrease, Ca 2+ dissociates from CaM before CaM dissociates from its binding partner [81]. Of course, the affinity of apo-CaM for CaMKII (1.45mM) should be insufficient to out-compete PP1-CaMKII binding (0.16μM). Yet because the 9-state-1-step model explicitly accounts for CaM-trapping by increasing the affinity of pThr-286 subunits for Ca 2+ /CaM by 1000-fold [73], the affinity of pThr-286 CaMKII for apo-CaM (1.45μM) is within one order of magnitude as that of PP1 (0.16μM). (In the original 2-state-2-step model, we provide for CaM-trapping implicitly through the 2-step CaM-binding.) Thus, apo-CaM may be able to compete with PP1 for CaMKII-binding, but it now remains for future experimental studies to directly quantify the kinetics and/or structure of the apo-CaM interaction with pThr-286 CaMKII. Determining the effect of CaM-trapping on sub-saturated Ca 2+ /CaM states is outside the scope of the current work. Notably, during dynamic Ca 2+ flux (S1 Appendix), we observe that non-apo-Ca 2+ /CaM states such as CaM 2C may contribute to maintenance of pThr-286. This indicates that even if structural studies reveal that apo-CaM is insufficient to out-compete PP1, it is possible that CaM-dependent PP1 exclusion could contribute at least somewhat to pThr-286 acquisition and maintenance during Ca 2+ stimulation. Clearly, further structural and kinetic studies of the CaM-CaMKII and PP1-CaMKII interaction are needed.
CaM-dependent PP exclusion could provide an added layer of robustness to similar mechanisms that may protect pThr-286 from de-phosphorylation. For example, Mullasseril et al. (2007) observe that endogenous, PSD-resident PP1 cannot de-phosphorylate CaMKII at pThr-286, whereas adding exogenous PP1 does cause de-phosphorylation [71]. The results by Mullasseril et al. suggest that endogenous PP1 is somehow sequestered by the PSD scaffold, and only upon saturation of this scaffold by exogenous PP1 does pThr-286 become de-phosphorylated. Our results indicate that perhaps in addition to saturating the PSD scaffold, the added exogenous PP1 could be out-competing CaM for binding to CaMKII, thereby terminating protection of pThr-286 by CaM. As another example, Urakubo et al. suggest that pThr-286 could be protected from PP activity by GluN2B binding, showing that GluN2B peptides are necessary for an apparent CaMKII bistability in vitro [72]. Notably, Urakubo et al., when using equal concentrations of CaM and CaMKII subunits, observe a slow and steady decline in pThr-286 upon kinase inhibition that has a similar rate as that seen in the Bradshaw study [80], taking place over the course of hours, which we only observe in our exclusive model. Urakubo et al. also mention that when basal levels of Ca 2+ are set below 0.2 μM, CaMKII autophosphorylation decreases to basal levels within 6 hours, perhaps indicating that GluN2Bdependent exclusion of PP1 is necessary but not sufficient for maintaining pThr-286 levels; sufficiently high levels of Ca 2+ /CaM may also be required. Overall, it seems scaffold-dependent sequestration of PP1 [71], GluN2B-dependent PP exclusion [72], and CaM-dependent PP exclusion could together provide considerable robustness of pThr-286 to phosphatase activity.

Simulation methods
In each MCell execution, proteins are instantiated at time zero having random positions within a 0.0328 μm 3 (0.0328 fL) cube, with each edge being 0.32μm in length. All proteins are described as three-dimensional volume molecules having the following concentrations: 1.52 μM CaMKII (30 holoenzymes, 18.24 μM and 360 individual CaMKII subunits), 22.8 μM CaM (450 discrete proteins), and 0.86 μM PP1 (17 discrete proteins). Because CaMKII particles are modeled using the specialized COMPLEX_MOLECULE syntax and MCell 3.3 does not accommodate diffusion for such particles, CaMKII is given no diffusion constant. In contrast, CaM and PP1 are simple volume-type molecules that move about the model space with a diffusion constant 6 × 10 −6 cm 2 /sec, so chosen to minimize the effects of any possible spatial localizations that may arise during a simulation. We emphasize that because this model does not explore spatial effects and, indeed, does not utilize a physiological spine geometry, using such a relatively fast diffusion parameter ensures that spatial effects do not confound our results. Future models exploring the spatial dependence of CaMKII holoenzyme activity in the dendritic spine may choose to adopt different diffusion parameters. All models are run at a time step of 0.1 μs for a total of either 20 or 120 seconds of simulation time, depending on the model variant.
MCell is a particle-based spatial-stochastic simulation engine. In a particle-based framework, individual protein species are modeled as discrete objects in space, rather than bulk/ well-mixed fluids. At each model timestep, MCell calculates each protein particle's subsequent diffusion distance and trajectory, in addition to the particle's probabilities for reacting with any nearby particles. More information about MCell's internal algorithms may be found at mcell.org and in publications such as those by Bartol et al. [51]. In short, the particle-based framework in MCell provides for spatial and stochastic considerations because each protein particle has unique spatial coordinates that proceed along random (stochastic) trajectories. Importantly, we assert that spatial-stochastic frameworks may be essential to characterizing CaMKII regulation, because 1) proteins in the spine are spatially organized and 2) protein copy numbers in the spine are low (tens to hundreds each), possibly invalidating the Law of Mass Action. Because MCell models are stochastic and change with each simulation, we average the output of many identical simulations. To ensure that the averaged output converges and is statistically significant, all model variants are repeated 50 times each.
CaM activation/inactivation is modeled by a pair of forcing functions which serve as a proxy for Ca 2+ flux. Both forcing functions are time-dependent square waves and inform the rates at which free CaM transitions between states (see Fig B in S1 Appendix) Eq 1 rapidly transitions all free CaM towards an active (Ca 2+ /CaM) state, and Eq 2 rapidly transitions all free CaM towards an inactive (apo-CaM) state. Time t iterates at 0.01sec intervals for the complete duration of a simulation. Eqs 1 and 2 therefore yield a peak width of 0.01sec regardless of frequency, which allows us to directly compare the effect of different Ca 2+ /CaM frequencies on CaMKII activity, without having to account for variable amounts of Ca 2+ /CaM exposure per pulse. In separate simulations without frequency dependence (i.e. Ca 2+ /CaM is continuously available to CaMKII), Eq 1 is adjusted to always fulfill the t = n i condition. Similarly, for pulse simulations in which Ca 2+ /CaM becomes withdrawn or blocked, Eqs 1 and 2 are given abbreviated time domains. All MCell code and associated files are available online at Github, the Purdue University Research Repository, and the University of Edinburgh Repository. Material will be made available upon publication.
Supporting information S1 Appendix. This document enumerates the model parameters, discusses combinatorial explosion, shows alternative visualizations of select data, and discusses the quantitative basis for PP1/CaM mutual exclusion from CaMKII binding. (DOCX)