Abstract
Activation of calciumcalmodulin dependent protein kinase II (CaMKII) during induction of longterm potentiation (LTP) is a series of complicated stochastic processes that are affected by noise. There are two main sources of noise affecting CaMKII activation within a dendritic spine. One is the noise associated with stochastic opening of Nmethyldaspartate (NMDA) receptor channels and the other is the noise associated with the stochastic reactiondiffusion kinetics leading to CaMKII activation. Many models have been developed to simulate CaMKII activation, but there is no fully stochastic model that studies the effect of noise on CaMKII activation. Here we construct a fully stochastic model to study these effects. Our results show that noise has important effects on CaMKII activation variability, with the effect from stochastic opening of NMDA receptor channels being 5–10 times more significant than that from stochastic reactions involving CaMKII. In addition, CaMKII activation levels and the variability of activation are greatly affected by small changes in NMDA receptor channel number at the synapse. One reason LTP induction protocols may require tetanic or repeated burst stimulation is that there is a need to overcome inherent variability to provide sufficiently large calcium signals through NMDA receptor channels; with meaningful physiological stimuli, noise may allow the calcium signal to exceed threshold for CaMKII activation when it might not do so otherwise.
INTRODUCTION
One of the most studied phenomena in neuroscience over the past three decades is longterm potentiation (LTP). LTP is a modification of synaptic strength originally discovered in the dentate and is the best model we have of how memory and learning occur (Bliss and Collingridge 1993; Malenka and Nicoll 1993). In the most studied form of LTP, high levels of calcium influx through Nmethyldaspartate (NMDA) receptor channels at synapses on dendritic spines activate calcium/calmodulindependent protein kinase II (CaMKII) to induce LTP (Lisman 2003; Lisman et al. 1997; Malenka and Bear 2004). The calcium that enters the spine head from the synaptic cleft binds to calmodulin to produce the calcium–calmodulin complex that can bind to individual subunits of CaMKII and activate them. When enough CaMKII subunits are activated, autophosphorylation can occur (Miller and Kennedy 1986), which traps calmodulin on a subunit and makes the subunit active long after the calcium signal is over (Meyer et al. 1992). The actions of activated CaMKII induce LTP (Fukunaga et al. 1993, 1996; Lisman et al. 2002).
There are many types of noise in the nervous system (Ermentrout et al. 2008; Faisal et al. 2008). Three main types of noise in neurons are thermal noise (Jaramillo and Wiesenfeld 1998; Johnson 1987; Shimozawa et al. 2003), ion channel noise (Clay and DeFelice 1983; White et al. 2000), and synaptic noise (Calvin and Stevens 1967; Fellous et al. 2003). Noise in neural systems has important effects on many neural processes, such as spike generation (Faisal and Laughlin 2007; Zeng and Jung 2004), information transmission (Bezrukov and Vodyanoy 1995; Hoch et al. 2003; Yu et al. 2001; Zeng et al. 2007), perception (Shimozawa et al. 2003), and the observed trialtotrial variability in LTP (Malenka and Nicoll 1993; Malinow and Tsien 1990). In particular, CaMKII activation, which leads to LTP induction, is a series of stochastic reactiondiffusion processes that are subject to noise. Others have shown that noise can affect the induction and stability of the CaMKII switch (Bhalla 2004a,b; Graupner et al. 2007; Hayer and Bhalla 2005; Miller and Wang 2006; Miller et al. 2005), which is implicated in LTP/longterm depression (LTD) and longterm memory (Lisman and Zhabotinsky 2001; Lisman et al. 2002). Thus it is important to understand the effect of noise on CaMKII activation.
There are many sources of noise that affect CaMKII activation within a dendritic spine. One source is the stochastic opening of NMDA receptor channels (Chen et al. 1999; Rosenmund et al. 1995). Another is the stochastic reactiondiffusion kinetics leading to CaMKII activation within a dendritic spine, which includes stochastic binding of calcium to calmodulin (Keller et al. 2008), stochastic binding of calmodulin to CaMKII (Byrne et al. 2009), stochastic calmodulin trapping, and stochastic calciumindependent autophosphorylation (Franks and Sejnowski 2004; Xia and Storm 2005). Models have been developed to simulate how the calcium signal leads to calmodulin trapping and CaMKII activation (Coomber 1998; Dosemeci and Albers 1996; Holmes 2000; Kubota and Bower 2001; Matsushita et al. 1995; Michelson and Schulman 1994; Okamoto and Ichikawa 2000; Zhabotinsky 2000; Zhabotinsky et al. 2006). In most of these models, the calcium signal was considered to be an arbitrary wave form, and binding of calcium to calmodulin, calmodulin trapping, and calciumindependent autophosphorylation were treated deterministically. However, because of the small numbers of some components in these reactions, these processes should be treated stochastically.
Here we use a stochastic model of reaction and diffusion at a synapse to determine NMDA receptor channel openings. These NMDA receptor channel openings are used in a neuron level model to compute calcium influx through NMDA receptor channels in dendritic spines following tetanic input to hundreds of synapses in a model of a fully reconstructed dentate granule cell. This calcium influx is used in a fully stochastic model of calcium dynamics in a dendritic spine to simulate binding of calcium to calmodulin, calmodulin trapping, and calciumindependent autophosphorylation. We find that CaMKII activation can be highly sensitive to stochastic openings of NMDA receptor channels. The effect from stochastic openings of NMDA receptor channels is almost 5–10 times larger than that from noise in the stochastic reactions leading to CaMKII activation. In addition, CaMKII activation is greatly affected by small changes in the number of NMDA receptors at a synapse.
METHODS
Generation of the calcium input signal
The calcium input signal is calculated with the use of synapse and neuron level models as described in Holmes (2000), except that here the opening of NMDA receptor channels is modeled stochastically at the synapse of interest. In the synapse level model, glutamate (2,000 molecules) is released into the synaptic cleft where it diffuses, gets pumped into neurons or glia by glutamate transporters, or binds to NMDA and AMPA receptors, leading to the opening of NMDA and AMPA receptor channels at the synapse. This process is computed deterministically at most synapses (Holmes 1995) but stochastically with MCell v. 2.50 (Stiles and Bartol 2001) at the synapse of interest. For a fixed number of AMPA and NMDA receptors, the locations of receptors and glutamate uptake transporters was identical across trials to eliminate variability of receptor location as a source of noise. The time course of the number of open NMDA and AMPA receptor channels is used at synapses in a neuron level model of a fully reconstructed dentate granule cell with voltagedependent conductances on the soma, axon, and proximal dendrites (Holmes and Levy 1997). Activation of 300 excitatory synapses on the dendritic spines in the middle third of the dendritic tree is modeled to represent strong medial perforant path (mpp) input. These synapses are activated with a single stimulus train of eight presynaptic spikes at frequencies of 10–400 Hz. The voltage computed in the neuron level model is used to determine the number of open NMDA channels (from the synapse level model) that are not blocked by magnesium. Blockage by magnesium is modeled deterministically except at the synapse of interest where it is modeled stochastically. Given a single channel conductance of 50 pS for NMDA receptor channels, the number of open, unblocked channels at the synapse, and the voltage at the synapse, the NMDA current at the synapse can be computed. The calcium portion of the NMDA current is computed with the constant field equation. This calcium portion of the NMDA current is used as the calcium signal in the spine model described below.
MCELL MODEL DETAILS
The parameters of the synapse model including the reaction schemes for NMDA and AMPA receptor channels were those used previously (see Holmes and Levy 1997 for reaction schemes and Holmes 1995 for cleft parameters). In MCell, the synaptic cleft was a box 8 × 8 × 0.02 μm with the vesicle release site located at the center of the top of the box and the postsynaptic density (PSD) located at the center of the bottom of the box. A synaptic vesicle 0.05 μm in diameter released 2,000 glutamate molecules into the synaptic cleft. The glutamate diffusion constant was 0.25 × 10^{−5} cm^{2}/s. The PSD radius was 0.1 μm. We chose NMDA and AMPA receptor densities to match different numbers of NMDA and 63 AMPA receptors. However, MCell 2 randomly assigns receptor number from the density values (and also location), and this means different seeds will produce different numbers of receptors for the same densities. We ran simulations with different seeds to find a set of seeds that would give us, for example, 14 and 63 NMDA and nonNMDA receptors, respectively. A seed that gave us the desired number of receptors was selected and was used for one time step to fix the locations of NMDA receptors, AMPA receptors, and glutamate uptake transporters. Making use of MCell's checkpoint feature, we chose a different seed for each of 30 simulation trials to simulate all subsequent time steps. We note that this did not fix the number of glutamate uptake transporters across simulations with different numbers of NMDA receptors because MCell assigns transporter number and location randomly. It was not possible to find seeds that gave the same exact uptake environment for every desired combination of NMDA and AMPA receptors. To make uptake characteristics as consistent as possible across all simulations, we modified uptake rate constant values. Specifically, we model neuronal and glial glutamate uptake with MichaelisMenten kinetics with fixed K_{m} and V_{max} values for each type of uptake as in Holmes (1995); this does not fix transporter density and transporter uptake rate independently, so in MCell, we fixed uptake transporter density and modified the two forward rate constants in the MichaelisMenten scheme to match V_{max} and K_{m} based on the number of transporter molecules MCell decided to choose for density for a particular random seed from the density value (backward rate constants were assumed to be zero).
COMPUTING THE CALCIUM COMPONENT OF THE NMDA CURRENT.
The number of open NMDA receptor channels as a function of time computed with MCell was read into the neuron level model for the synapse of interest. When a channel became open, it was assigned randomly to be either open or blocked according to the probabilities computed from kinetic rates reported by Ascher and Nowak (1988). At each time step, each open NMDA channel and each blocked NMDA channel was interrogated to determine whether the channel would change state, given the voltage at the synapse and the probabilities of switching states. If the MCell input indicated that a channel closes, probabilities based on the rate constants were again used to decide if a blocked channel closed or an open channel closed. Extracellular magnesium concentration was 1.2 mM. The dentate granule cell model computes I_{NMDA} = g_{NMDA}(V – V_{NMDA}), where g_{NMDA} is the number of open, unblocked NMDA receptor channels multiplied by 50 pS, the single channel conductance. Given the ion concentrations, V_{NMDA} was +1.3827 mV, slightly larger than the value of 0 mV usually used. To separate out the calcium component, we use the procedure from Holmes and Levy (1990), which is based on the analysis by Mayer and Westbrook (1987) using the GoldmanHodgkinKatz (GHK) constant field current equations and ion activities. From Eqs. 3 and 4 from Mayer and Westbrook (1987), we get
Mayer and Westbrook (1987) report that P_{Na}:P_{K}:P_{Ca} = 1:1:10.6. Substituting in the preceding equations, we get I_{NMDA} = P_{Na}(c + d + 10.6e). This is solved for P_{Na}, and I_{Ca} = 10.6P_{Na}e. Mayer and Westbrook (1987) found better fits to the data when they used ion activities instead of concentrations. Thus the concentrations given in the preceding equations were multiplied by their activity coefficients in our computations. We use 10 mM Na^{+}, 140 mM K^{+}, and 20 nM Ca^{2+} inside the cell and 145 mM Na^{+}, 2.5 mM K^{+}, and 2 mM Ca^{2+} outside. Activity coefficients are 0.75 for Na^{+}, 0.75 for K^{+}, and 0.19 for Ca^{2+}, and temperature is 37°C making RT/F 26.73.
Stochastic simulation of reaction and diffusion
Stochastic reactions in a homogeneous volume are often modeled with the Gillespie algorithm (Gillespie 1976, 1977). However, the living cell is usually a spatially complicated system that is not homogeneous. Therefore the Gillespie method has to be adjusted to fit the heterogeneous system. A typical solution for spatially heterogeneous coupled reactiondiffusion systems is to divide the system into many voxels, which can be considered as homogeneous (Gillespie 1976). The homogeneous Gillespie's method can be applied to each individual voxel. We will briefly describe the Gillespie algorithm in one volume and show how we extend it to heterogeneous volumes.
First consider reactions within one homogeneous volume. A system consisting of coupled chemical reactions may contain reactions (R_{1},_{,} …, R_{K}), each having the following stoichiometric form
If there are X_{1}, …, X_{N} molecules of species S_{1}, …, S_{N}, the number of distinct reactant combinations, h_{μ}, for reaction μ is
Given the number of distinct reactant combinations and the stochastic rate parameter c_{μ}, the probability of reaction μ occurring in unit time is a_{μ} = h_{μ}c_{μ}. The sum of the probabilities over all reactions is
According to the Gillespie algorithm, one can determine the time of the next reaction and the identity of the next reaction by drawing two uniformly distributed random numbers r_{1} and r_{2}. The time of the next reaction τ is
The simulation is advanced a time step τ and the numbers of molecules X_{1}, …, X_{N} of species S_{1}, …, S_{N} are updated according to the stoichiometry of reaction μ.
Although this algorithm is straightforward to implement, it is valid only for a homogeneous wellmixed volume. To use this algorithm for an inhomogeneous volume, the volume is divided up into voxels that individually are homogeneous and well mixed. Gillespie suggested a method based on Fick's first law of diffusion j = −D▿C, where flux j can be interpreted as the number of diffusing ligands that cross the unit surface in unit time. The diffusion of a molecule of species α from voxel V_{l} into its neighbor V_{l+1} can be simulated as a special reaction consisting of two simultaneous subreactions: S_{α} → * (destroyed) in V_{l} and * → S_{α} (created) in V_{l+1}. By Fick's first law of diffusion, the probability coefficient d_{α} can be obtained as
The a defined by Eq. 4 is modified to include the probabilities for all reactions and all diffusions in all voxels
Hence, we can outline the stochastic numerical simulation algorithm.

1) Initialization: initialize all relative variables, including calculating the “reaction rates” for all “reactions” in all voxels, and the initial numbers of molecules in all voxels. Calculate the variable h_{μ}, d_{α} in each voxel and calculate a as defined above. Also set initial time t = 0.

2) Advance the simulation and update the system states: first, generate one pair of random numbers (r_{1}, r_{2}), calculate τ, and determine the (ijk)_{th} voxel where the “reaction” R_{μ} occurs; second, advance t by τ and update the numbers of molecules of those species involved in “reaction” R_{μ} in the (ijk)_{th} voxel (note: if it is a diffusion, the quantity of the corresponding molecules in the appropriate neighboring voxel, where the molecules diffuse also has to be updated). Finally, recalculate the variable a based on the new quantities of molecules in the relevant voxels.

3) If t does not reach the ending point, write out t and the molecule quantities of interest and return to step 2; if t reaches the ending point or all molecules have vanished, terminate the simulation.
The method only requires two random numbers for each time step; hence the computation for random number generation encountered in traditional simulations is reduced dramatically. This allows the simulation to include a large number of subvolumes (voxels) for more realistic models.
To verify that the model was working properly, results were compared with those obtained with a mixed stochasticdeterministic model described earlier (Holmes 2000).
Dendritic spine model
The dendritic spine model used here is extended from that used previously (Holmes 2000). Taking into account calcium diffusion, pumping, and binding to buffers, this model calculates spine head calcium concentration after calcium influx through NMDA receptor channels. The spine head is chosen to be of longthin type. The spine head is a 0.5 × 0.5 × 0.55μm rectangular box, and the spine neck is a 0.1 × 0.1 × 0.8μm rectangular box, connecting to the underlying dendrite, which is represented by a 2 × 1 × 1μm rectangular box. Because CaMKII and calcineurin molecules are mostly localized within the outer 50 nm of the spine head, the outermost layer of the spine head is filled with 50 × 50 × 50nm voxels, and the rest of the spine head is filled with 0.1 × 0.1 × 0.1μm voxels. There are a total of 225 voxels in the spine head. The spine neck is divided into 8 0.1 × 0.1 × 0.1μm voxels, and the underlying dendrite is divided into 16 0.5 × 0.5 × 0.5μm voxels.
In the spine, calcium influx through NMDA receptor channels enters the spine head uniformly through the middle 16 voxels in the outer layer of the spine head, representing the region of the postsynaptic density. Calcium can bind to calbindin or calmodulin, and calmodulin with zero to four calcium ions bound can bind to CaMKII, calcineurin, or neurogranin. Calcium and calmodulin can diffuse freely, but CaMKII and calcineurin are restricted to the outermost layer of the spine head. There are 83 CaMKII holoenzymes (each with 2 rings of 6 subunits each) and 200 calcineurin molecules whose primary function in these short simulations is to compete with CaMKII for calmodulin. Neurogranin is present in all areas of the spine and dendrite with a total concentration of 26.7 μM. Total calbindin concentration is 40 μM, and total calmodulin concentration is 13.3 μM. Recently CaMKII was shown to be activated by calmodulin having only two bound Ca^{2+} ions (Shifman et al. 2006); therefore in the model, a CaMKII subunit bound with CaMCa_{x}, where x is ≥2 is considered to be in the “bound” state. A bound subunit can become autophosphorylated and reach the “trapped” state (CaMCa_{x} bound and phosphorylated) only if its immediate neighbor on the same ring of the CaMKII holoenzyme is also in the bound or trapped state (or with lower probability, in the autonomous or capped states). If a trapped subunit loses CaMCa_{x}, it is considered to be in the autonomous state and, from there, autophosphorylation can occur in the calmodulin binding site to bring the subunit into the “capped” state if the neighboring subunit in the holoenzyme ring to the right or left is bound, trapped, autonomous, or capped.
RATE CONSTANT DETAILS
The two calcium buffers in the model are calbindin and calmodulin. The concentration of calbindin in adult dentate granule cells has been estimated to be 40 μM (Müller et al. 2005) and that is the concentration used in the model. Calbindin can bind four Ca^{2+} ions; two sites are high affinity sites and two are low affinity. On and off rates in the model were 5.5 μM^{−1} s^{−1} and 2.6 s^{−1} for the highaffinity site and 43.5 μM^{−1} s^{−1} and 35.8 s^{−1} for the lowaffinity site (Nagerl et al. 2000; with modification of on rates suggested by Berggard et al. 2002). Calmodulin also binds four Ca^{2+} ions, two on the Nlobe and two on the Clobe. The K_{d} values used were 10 and 1.2 μM for the Clobe and 50 and 10 μM for the Nlobe (Linse et al. 1991). The on and off rates are shown in the back frames of Fig. 1, A–C. Dissociation rates are based on data from Gaertner et al. (2004) and, together with the K_{d} values, these determined the on rates. Diffusion coefficients were 0.223 for calcium (Allbritton et al. 1992), 0.0223 for calmodulin (Schmidt et al. 2007), and 0.02 for calbindin (Schmidt et al. 2005) (all in units of μm^{2}/ms). Most of the calmodulin is bound to neurogranin at rest.
Rate constants for Ca^{2+} binding to CaMKIICaMCa_{x} and CaMCa_{x} binding to CaMKII are given in Fig. 1A. For Ca^{2+} binding to CaMKIICaMCa_{x}, the dissociation rates are from Gaertner et al. (2004); the K_{d} values were assumed to be 10fold higher affinity than calcium binding to CaMCa_{x} without the kinase, and the on rates were computed from the K_{d} and the off rates. These rates are shown in the front frame of Fig. 1A. For CaMCa_{x} binding to CaMKII, the K_{d} and dissociation rate for CaMCa_{4} binding to CaMKII are based on values from Meyer et al. (1992). The other values on the slanted lines in the figure are computed from thermodynamic considerations, given the remaining rate constants.
Rate constants for Ca^{2+} binding to calcineurin bound with CaMCa_{x} and CaMCa_{x} binding to calcineurin (CaN) are given in Fig. 1B. For Ca^{2+} binding to CaNCaMCa_{x}, the K_{d} values of 0.1, 0.06, 0.3, and 0.08 μM were taken from Stemmer and Klee (1994). The on and off rates were assumed to be in the same range as those with CaMKII. These rates are shown in the front frame of Fig. 1B. For CaMCa_{x} binding to CaN, the K_{d} and on and off rates for CaMCa_{4} binding to CaN were taken from Quintana et al. (2005). The other values on the slanted lines in this figure are computed from thermodynamic considerations, given the other rate constants.
Rate constants for Ca^{2+} binding to neurogranin bound with CaMCa_{x} and CaMCa_{x} binding to neurogranin (Ng) are given in Fig. 1C. For Ca^{2+} binding to NgCaMCa_{x}, the dissociation rates were taken from Gaertner et al. (2004) and the on rates were assumed to be the same as calcium binding to calmodulin without neurogranin. These rates are shown in the front frame of Fig. 1C. For CaMCa_{x} binding to Ng, we assumed that CaM binding to Ng had a K_{d} of 3 μM (reports range 1–5 μM; see Zhabotinsky et al. 2006). The on and off rates were chosen from thermodynamic considerations but were similar to those in Kubota et al. (2007).
CaMKII subunits in the bound state as defined above can become trapped or phosphorylated if a neighboring subunit on the same ring of the holoenzyme is also bound, trapped, and with 0.4 times lower probability, autonomous or capped. The trapping rate constant is 10 s^{−1}, which is much higher than in our previous models but seems more consistent with data from Bradshaw et al. (2003). We do not model dephosphorylation in detail here; the rate used is 0.003 s^{−1}, making this a rare event within the limited time course of the results shown. Subunits in the trapped state can gain or lose Ca^{2+} ions. Our rate constants for this are given in the back frame of Fig. 1D. The dissociation rates were taken from Gaertner et al. (2004), and the association rates were chosen to make the K_{d} of each step the same as with the unphosphorylated kinase. Rate constants for the transition from trapped states to the autonomous state are also shown. Rates were chosen to be consistent with Meyer et al. (1992) (Fig. 3), although these reactions were not common within the time course of the current simulations.
RESULTS
Combined effect of two kinds of noise
We first studied the combined effect of noise from both stochastic opening of NMDA receptor channels and stochastic CaMKII activation. Opening of NMDA receptor channels at the synapse of interest was modeled stochastically with MCell v. 2.50 (Stiles and Bartol 2001), and CaMKII activation within the dendritic spine was also simulated stochastically. Opening of NMDA receptor channels on other synapses was modeled deterministically. An eightpulse tetanus at 100 Hz was applied to 300 synapses in a model of a dentate granule cell to determine the voltage change at the synapses. At each synapse, each pulse in the tetanus was modeled as the release of 2,000 glutamate molecules. Each deterministic synapse had 63 AMPA and 14 NMDA receptor channels. The stochastic synapse had 63 AMPA receptors and 10–26 NMDA receptors depending on the simulation trial group. We used MCell results together with the voltage from the neuron level model to generate 30 trials of calcium signals, which were used in the spine model to compute levels of CaMKII activation. In the 30 trials, all of the simulation conditions were the same, including the locations of the NMDA receptors, AMPA receptors, and glutamate uptake transporters, except that different random seeds were used in all MCell and stochastic CaMKII activation simulations.
The number of NMDA receptor channel openings computed stochastically with MCell after an eightpulse, 100Hz tetanus showed considerable variability across trials as shown in Fig. 2. Trials with maximal opening and minimal opening, representing extreme cases encountered among 30 trials when there were 19 NMDA receptors at the synapse, are compared in Fig. 2A. In Fig. 2A (top), the peak number of open NMDA channels was eight, and there were long stretches when there were always one or two open channels. In Fig. 2A (middle), the peak number of open channels was three, and stretches where there was at least one open channel were comparatively short. NMDA receptor channels that reach the open state are subject to voltagedependent magnesium block. Taking into account the voltage computed in the neuron level model, the peak number of open and unblocked channels was three and two for these two trials, respectively, as shown in Fig. 2B (top and middle). The bottom traces in Fig. 2 show the average number of open NMDA channels both with and without magnesium block. These traces show some late channel openings commonly seen in these trials. We sometimes saw isolated openings 900–1,000 ms after the tetanus.
The variability in the number of open and unblocked NMDA receptor channels translated into very different calcium signals in the outer layer of the dendritic spine head as shown in Fig. 3A. Here we show the calcium concentration for the two extreme cases from Fig. 2 and, for comparison, the calcium concentration averaged over all 30 trials. The large difference in calcium concentration was expected because total calcium influx varied over a two to threefold range among the pictured trials (11,000–26,000 Ca^{2+} ions). This difference in the calcium signal caused a huge difference in the number of activated CaMKII subunits, as shown in Fig. 3B. Results from all 30 trials, as well as the average, are shown in this figure, and they show that the two extreme cases from Fig. 2 are not isolated examples but reflect the inherent variability in CaMKII activation.
Because the induction of LTP is known to require a highfrequency tetanus, we next sought to examine the variability in CaMKII activation for different tetanus frequencies. We produced plots of CaMKII activation similar to Fig. 3B for 10 to 400Hz stimulation, which are summarized in Fig. 4A. These results again show that the effect of noise on CaMKII activation can be substantial. Perhaps surprisingly, a large number of activated CaMKII subunits occurred for two 20Hz trials, and some small numbers of activated subunits occurred in some 100Hz trials. The variability among the 20Hz trials and their time courses are pictured in Fig. 4B. Nevertheless, the mean and median number of activated CaMKII subunits generally increased with frequency in the expected manner for 10–100 Hz. The illustrated example shows a dip for 200 and 400Hz tetani trial groups, but this was not a consistent finding; results with 18 NMDA receptors at the synapse showed constant means and medians for 100 to 400Hz trial groups.
NMDA receptor number provides an additional source of variability
To this point, the illustrated simulations always used 19 NMDA receptors at the stochastic synapse. How would small changes in the number of NMDA receptors at the synapse affect CaMKII activation? To answer this, we used eightpulse 20 and 100Hz tetani with different numbers of NMDA receptor channels to generate the input calcium signals, and then we computed CaMKII activation. We simulated 30 trials each for 10–26 NMDA receptor channels at each tetanus frequency.
We found that, although the general trend, as expected, was that more CaMKII subunits were activated as the number of NMDA receptor channels or the tetanus frequency was increased, there was considerable variability among runs as shown in Fig. 5, A and B. Means and medians only are plotted in Fig. 5C to show trends more clearly. In the 100Hz tetanus results, the mean number of activated CaMKII subunits increased fourfold and the median increased threefold as the number of NMDA receptors was raised from 14 to 18; the mean increased another twofold and the median threefold as the number of NMDA receptors was raised from 18 to 24. In the 20Hz tetanus results, increases in the number of activated CaMKII subunits with increasing NMDA receptor number were much more gradual; nevertheless, there were a number of trials at this tetanus frequency at nearly every NMDA receptor number where CaMKII activation exceeded the corresponding mean or median at 100 Hz.
Individual effect of two kinds of noises
How much of the observed variability comes from stochastic opening of NMDA receptor channels and how much comes from stochastic CaMKII activation reactions? We computed the calcium influx generated by an eightpluse 100Hz tetanus at a synapse with 19 NMDA receptor channels and used this same calcium signal in the dendritic spine model but with different random seeds for the CaMKII activation reactions. The calcium signal chosen was the trial that had calcium influx closest to the mean of 30 trials. The results suggest that CaMKII reactions provide a much smaller source of variability than NMDA receptor channel openings (Fig. 6).
We next sought to quantify this result and extend it to include the effect of different numbers of NMDA receptors at the synapse. To quantify the combined effect of both types of noise on CaMKII activation variability, we first measured the relative fluctuation √<̅(̅N̅ ̅−̅ ̅<̅N̅>̅)̅2̅>̅/̅<̅N̅>̅ of the total number of activated CaMKII subunits, where N is the sum of the number of bound, trapped, autonomous, and capped subunits in each trial at 1 s and <N> is the average over 30 trials. To separate the two sources of noise, we compared results computed with stochastic models of NMDA receptor channel openings and a stochastic model of CaMKII activation to results computed, as for Fig. 6, with the trial whose NMDA receptor channel openings and subsequent calcium influx was closest to the mean of 30 trials and stochastic models of CaMKII activation. This was done for NMDA receptor numbers of 12–26. Results are shown in Fig. 7.
We see from Fig. 7 that the relative fluctuation of the number of activated CaMKII subunits attributed to stochastic CaMKII activation alone is just 12–28% of that coming from both stochastic opening of NMDA receptor channels and stochastic CaMKII activation, with the percentages being 12–15% for larger numbers of NMDA receptors. This was to be expected because, in the spine head, there are more CaMKII subunits (∼1,000 in these simulations) than NMDA receptor channels (∼20). As the number of NMDA receptor channels was increased, the relative fluctuation in the number of activated CaMKII subunits decreased. For example, when NMDA receptor channel number was 12, the relative fluctuation from stochastic opening of NMDA receptor channels and stochastic CaMKII reactions was 110%; however, when receptor number was 26, the relative fluctuation was reduced to 60%. Similarly, the relative fluctuation when only CaMKII reactions were modeled stochastically dropped from 22 to 8% over the same span of NMDA receptor number. One reason for this phenomenon is that increasing the number of NMDA receptor channels will decrease the fluctuation of the number of open NMDA receptor channels, which decreases the fluctuation of calcium influx, resulting in a decrease in the fluctuation in the number of activated CaMKII subunits. We note that the total number of activated CaMKII subunits in these simulations was low compared with the number of available subunits, especially when the number of NMDA receptors was small. This contributed to the high relative fluctuation shown in Fig. 7. As the number of NMDA receptor channels was increased and the number of activated CaMKII subunits became significant relative to the total number of subunits, the fluctuation naturally decreased.
DISCUSSION
It is well accepted that CaMKII activation is necessary and sufficient to induce some forms of LTP (Hayashi et al. 2000; Lledo et al. 1995). CaMKII activation occurs as a series of complicated stochastic processes, and there are a number of sources contributing to CaMKII activation variability within a dendritic spine (Chen et al. 1999; Franks and Sejnowski 2004; Rosenmund et al. 1995; Santucci and Raghavachari 2008; Xia and Storm 2005). Here we combine a fully stochastic synapse level model with a neuron level model and a fully stochastic model of reactiondiffusion within a dendritic spine to focus on two sources of noise: stochastic opening of NMDA receptor channels and stochastic CaMKII activation within a dendritic spine. Our results show that noise has important effects on CaMKII activation variability. Of the two sources of noise examined here, noise involved with stochastic opening of NMDA receptor channels was 5–10 times more important than noise in stochastic CaMKII activation reactions. The reason for this is that there are many more CaMKII subunits than NMDA receptor channels. Here we had ∼1,000 CaMKII subunits (83 holoenzymes) in one dendritic spine compared with ∼10–20 NMDA receptor channels (Racca et al. 2000). According to the square root law of statistical theory, the precision of the sample mean improves with the square root of the sample size. Applying the numbers given above, we would expect that the variability in CaMKII activation should be the square root of 1000/20 or 7 times larger from the stochastic opening of NMDA receptor channels than from stochastic CaMKII activation kinetics. This is exactly the relationship we obtained in the simulations when NMDA receptor number was 18 or larger. With smaller NMDA receptor numbers, the ratio we found fell to 3–5, mostly likely the result of having very few CaMKII subunits becoming activated in these cases. We want to note that, in early simulations in which we excluded calbindin, a larger fraction of the CaMKII subunits became activated for all numbers of NMDA receptors as would be expected; qualitative results were similar to those shown, and the expected sevenfold relationship was also found. If the number of holoenzymes was much less than that used in the model, say 30, then following the above formula, we would expect variability in CaMKII activation to be only 4 times larger from stochastic opening of NMDA receptor channels than from stochastic CaMKII activation kinetics. In this case, fluctuation because of CaMKII activation kinetics would become more important, but fluctuation because of stochastic NMDA channel openings would still dominate.
There are many sources of noise that we have not explored in this study. These can be divided into 1) additional factors to be considered and 2) variations in parameter values chosen. Additional factors to be considered include variations or heterogeneity in receptor type (i.e., NR1/2A vs. NR1/2B) (Chen et al. 1999; Santucci and Raghavachari 2008), nonuniform probability of vesicle release (Hanse and Gustafsson 2001; Rosenmund et al. 1993), stochasticity in other CaMKII activation reactions, particularly gating of phosphatase activation/deactivation (Bhalla 2002a,b); variation in spine shape (Holmes 1990; Schmidt and Eilers 2009), variation in the site of vesicle release (Uteshev and Pennefather 1997), variation in the location of NMDA receptors in the PSD (Santucci and Raghavachari 2008), and variation with the use of different proposed NMDA kinetic reaction schemes (Banke and Traynelis 2003; Popescu and Auerbach 2003; Yuan et al. 2005). As for parameter values, possible sources of variability include: the number of CaMKII holoenzymes in the spine head (Peterson et al. 2003), the diffusion rate of glutamate in the cleft (Holmes 1995; Rusakov and Kullmann 1998; Ventriglia and Maio 2000), the diffusion rate of calcium inside the spine, glutamate uptake rates, calcium pump rates, number of glutamate molecules in a vesicle (Burger et al. 1989; Holmes 1995; Riveros et al. 1986) and synaptic cleft width (Rusakov and Kullmann 1998; Ventriglia and Maio 2000). We suspect that noise from these sources will be less than that from stochastic variations in NMDA receptor channel openings, but verification of this awaits future work.
What are the consequences of noise for LTP induction and learning and memory? Because of noise in NMDA receptor channel openings, it may be necessary to repeat a tetanus or have repeated short bursts of stimulation to reliably kick CaMKII activation to a high level (Li and Holmes 2000). This may be particularly important during development when the number of NMDA receptors at a synapse may be small, making variability of CaMKII activation large. Only those synapses that receive repeated stimulation or certain patterned stimulation can reliably overcome the effects of noise. In addition, calcium influx is very sensitive to the number of NMDA receptors at the synapse and the peak number of open NMDA receptor channels is small (Dalby and Mody 2003; Nimchinsky et al. 2004). Although the number of NMDA receptor channels at synapses has been found to be far less variable than the number of AMPA receptor channels (Takumi et al. 1999), the results here show that even small differences in the number of NMDA receptor channels per synapse can have large effects on CaMKII activation. With physiological patterns of stimulation, noise may allow the calcium signal to exceed threshold for CaMKII activation when it would not do so otherwise. Noise may be an essential feature for selecting synapses to undergo potentiation. It is worth noting that a noisedriven stochastic selection mechanism for potentiating synapses has been found optimal for storing multiple memories in network models (Fusi 2002).
GRANTS
This work was supported by National Institute on Alcohol Abuse and Alcoholism Grant AA14294 via the Collaborative Research in Computational Neuroscience program and by the Scientific Research Fund of the Hunan Provincial Education Department (07B075).
ACKNOWLEDGMENTS
We thank L. M. Grover and J. AmbrosIngerson for comments and suggestions on drafts of this manuscript. We acknowledge the work of Y. Li during early stages of this project through 2003 and also thank N. Desmond for the morphology of the dentate granule cell used in the simulations.
 Copyright © 2010 The American Physiological Society