Structural plasticity controlled by calcium based correlation detection

The connectivity structure of the cortex was found to be surprisingly dynamic in vitro and in vivo (Bonhoeffer and Yuste, 2002 ). Synapse formation and elimination exhibit a marked dependence on spiking activity, where higher activity promotes synapse formation (Le Be and Markram, 2006 ). Based on geometric considerations, Stepanyants et al. (2002) , Chklovskii et al. (2004) and Stepanyants et al. (2007) suggest as a basic design principle of the cortex the potential of any pair of neurons to form a connection on small length scales of a few hundred micrometers together with an activity dependent selection mechanism. This idea is supported by direct observation of spines approaching presynaptic partners in a promiscuous manner (Kalisman et al., 2005 ). Of these structurally possible (potential) synapses, only a small fraction (0.12–0.34) is actually realized (Stepanyants et al., 2002 ), and transitions from potential to actual synapses are observed in vitro at rates of up to 1.2 per cent per hour during increased spiking activity (Le Be and Markram, 2006 ). These newly formed, immature synapses are weaker than mature ones. Synapse pruning mostly affects weak, but already mature synapses. The relation of synaptic strength to synapse formation and synaptic death indicates that long term potentiation (LTP) and synapse formation may be controlled by similar mechanisms and the same may hold for long term depression (LTD) and synapse pruning. The observation of increased connectivity after prolonged spiking activity also prompts for a synapse pruning mechanism, which restores the typical synapse density. In support of this hypothesis, not only in the developing cortex, but also in adult monkeys in vivo, synaptic boutons emerge and disappear with rates of 7 per cent per week (Stettler et al., 2006 ).

The functional role of synapse formation and elimination is not well understood, but theoretical studies have exposed the benefits of the capability of structural remodeling (Chklovskii et al., 2004 ; Stepanyants et al., 2002 ). Due to the fact that the number of potential presynaptic partners of a neuron easily exceeds the actual number at any point in time by an order of magnitude, each synapse carries three to four bits of structural information. In the framework of spiking associative memories this beneficial functional role of structural plasticity can be exploited to increase the memory capacity (Knoblauch, 2006 ; Knoblauch et al., 2007 ). Therefore, synaptic formation and elimination is a candidate process for long term information storage in cortical networks during development and adulthood alike.

The microscopic mechanisms leading to the formation of new synapses and to the elimination of existing ones have not yet been completely revealed. However, there is some evidence, that newly formed synapses are created in an intermediate, silent state (Cohen-Cory, 2002 ; Kalisman et al., 2005 ). These frequently encountered silent synapses lack AMPA receptors but have NMDA receptors (Atwood and Wojtowicz, 2004 ). They bear a high potential for remodeling the neural circuit, since they can easily be converted into active synapses. The most probable mechanism is the translocation of AMPA receptors into the postsynaptic density (PSD). This is also observed (Cohen-Cory, 2002 ; Shi et al., 1999 ) during LTP for which NMDA receptor activation is a necessary precondition.

There is strong indication that synapse formation and synapse pruning are controlled in a similar way as LTP and LTD (Lüscher et al., 2000 ). Furthermore, Le Be and Markram (2006) found that the same antagonists prevent LTP and synapse formation, and conclude that the underlying mechanisms may be similar. Calcium entering the postsynaptic site through NMDA receptors is a probable messenger causing synapse maturation (Atwood and Wojtowicz, 2004 ). For the induction of LTP and LTD, which in parts come about by modulating the number of postsynaptic AMPA receptors, Cormier et al. (2001) showed that both depend in a threshold like fashion on the intracellular calcium peak amplitude: higher amplitudes lead to LTP and lower ones to LTD.

Spike pairing experiments showed the calcium signal in spines to depend on the correlated activation of the pre- and postsynaptic site: Nevian and Sakmann (2004 , 2006 ) found a marked dependency on the relative timing and independence on the position of the spine along the dendrite, making the NMDA receptor mediated calcium level in a spine an attractive candidate substrate to convey information about the correlation between presynaptic and postsynaptic activity.

For the calcium signal in a spine to control synaptic maturation or synaptic death it must activate a downstream signaling cascade. Calcium/calmodulin dependent kinase II (CaMKII) is a calcium activated kinase, which is crucial for the LTP of a synapse. In its activated state it can phosphorylate several structures among them AMPA receptors which respond with increased conductivity. There is also evidence, that CaMKII is involved in the insertion process of new AMPA receptors into the PSD. This causes LTP or turns a silent synapse (only having NMDA receptors) into an active one having AMPA and NMDA receptors. There is recent evidence from detailed biophysical modeling studies (Graupner and Brunel, 2007 ) that the activity dependent calcium influx can activate CaMKII in a bistable fashion and hereby explain spike timing dependent synaptic plasticity (STDP). For a recent review of phenomenological models of STDP see Morrison et al. (2008) .

CaMKII forms holoenzymes of two ring molecules consisting of six subunits each. A subunit can either be active or inactive. Transitions between the inactive and the active state are triggered by calcium signals of different amplitudes. Short and weak calcium signals typically lead to activation of a single subunit by binding calcium or calmodulin to it. After the calcium level has dropped, unbinding of calcium and hence deactivation occurs within 0.1–0.2 s. At larger calcium concentrations an active subunit of the molecule (to which calcium is already bound) can phosphorylate the neighboring subunit. This only requires one additional calcium molecule to bind to the second subunit to expose its phosphorylation site. So phosphorylation can propagate along the ring and the molecule remains active even after calcium has returned to the resting level. At resting calcium concentrations, protein phosphatase 1 (PP1) can dephosphorylate an active subunit, but a neighboring active site can immediately rephosphorylate it again. This regenerating effect explains the long time scales of several minutes for the deactivation of groups of active CaMKII molecules. Even longer time scales of hours of persistent activity are found at resting calcium concentrations in the special chemical environment of the PSD, where the concentration of PP1 is low compared to the number of CaMKII subunits. For a comprehensive review see Lisman et al. (2002) . Detailed biophysical simulations (Miller et al., 2005 ) have confirmed bistability between an active and an inactive state of whole populations of approximately 20 holoenzymes. This effect is due to saturation of the phosphatase in the active state. The study found life times of both states in the range of 100 years. The time scale of the attractor dynamics is on the order of tens of minutes, but for strong fluctuations of the calcium signal the bistability vanishes. Calcium can not only activate CaMKII (either directly or via calmodulin), but also protein phosphatases like calcineurin and protein phosphatase 1, which dephosphorylate CaMKII. These phosphatases have a higher affinity to calcium than CaMKII. Therefore, they become active at lower calcium concentrations and counteract the phosphorylation of CaMKII. This is in agreement with the finding that LTD in CA1 dendrites can be induced if the calcium concentration is below 180 nM, while LTP requires it to exceed 540 nM (Cormier et al., 2001 ). For a review see Cavazzini et al. (2005) .

To our knowledge, previous models for structural plasticity either used simplified neuron models (Butz et al., 2008 ; Dammasch et al., 1986 ) or plasticity rules depending on the firing rate alone and not taking into account formation and death of individual synapses (van Ooyen et al., 1995 ). Consequently correlation dependent structure formation is outside the scope of these works. In this modeling study, we investigate how the biologically known pathways outlined above interplay to achieve a mechanism capable of detecting correlation between the presynaptic and the postsynaptic spiking activity. We focus on the calcium control hypothesis: the calcium signal mediated by NMDA receptors is the beginning of the signaling cascade. The main features of NMDA receptors entering our model are: (1) Their fast binding to glutamate followed by slow unbinding. (2) The quasi-instantaneous removal of the magnesium block upon postsynaptic depolarization to open the channel. In these assumptions, our model is similar to previous work by Shouval et al. (2002) on a mathematical model to explain spike timing dependent plasticity (STDP) based on the properties of NMDA receptors and the calcium control hypothesis. In contrast to their work, we assume the postsynaptic depolarization by the backpropagating action potential (bpAP) to be a short event.

The next stage of the signaling pathway in our model is a calcium controlled bistable effector molecule, like e.g. CaMKII. The important properties for our model are: (1) The long time constants of sustained activation of each individual molecule by high calcium concentrations (bistability). We are interested in the regime, where calcium fluctuations dominate the activation dynamics and the slower attractor dynamics causing bistability of the whole population of molecules is negligible. (2) The ability of the kinase to influence synaptic plasticity via AMPA receptor insertion. We assume a minimum amount of the kinase to be necessary for promoting synapses from silent to functional and we assume that there a minimal amount of active kinase is required to prevent synapse death. (3) Two disjoint ranges of calcium concentration that control the transitions between the inactive and the active state of each molecule. (4) The relatively low number of molecules, which makes a statistical description essential. We derive a model for correlation detection based on this pathway and investigate its dynamics under realistic conditions of irregular spike trains. We show that the model represents a viable mechanism to sense correlation between the pre- and postsynaptic activity. Controlling synaptic pruning, it can implement a firing rate homeostasis. Cohen-Cory (2002) already suspected, that activity-dependent remodeling selectively stabilizes coactive incoming synapses and destabilizes others. We demonstrate that such cooperation and competition between synapses naturally emerges from the microscopic model and that a neuron can learn the correlations between neighboring inputs.

In the section “Spike Time Dependence of Postsynaptic Calcium Concentration” we explain the origin of spike timing dependence of the postsynaptic calcium concentration, mention the main findings of recent imaging experiments, and develop a model for the peak amplitude of the postsynaptic calcium signal. We show it to qualitatively reproduce the experimental findings. In the section “Ca²⁺ Transients Caused by Correlated Irregular Spiking” we show that for irregular spiking activity this model predicts a distinct dependency of the observed postsynaptic calcium signal on the correlation between the presynaptic and the postsynaptic spiking events. The section titled “A Counter for Correlated Events” derives a biologically motivated model of a mechanism to “count” correlated events and therefore to assess the degree of correlation between the presynaptic and the postsynaptic activity. The section titled “Rate Homeostasis by Synaptic Pruning” shows that controlling synaptic pruning by this correlation measure can act as a rate regulation for the postsynaptic neuron at low rates. In the section “Cooperation and Competition by Spatial Input Correlation” we demonstrate that cooperation and competition between synapses depends on the correlation between neighboring inputs and that a synaptic pruning process manifests these input correlations in the resulting network structure. The last section discusses our results.

All simulations were carried out with the NEST simulation software (Gewaltig and Diesmann, 2007 ) using the computationally efficient implementation of synaptic maturation and death provided in the section “Algorithmic Implementation of Synapse Maturation and Synapse Death” in Appendix. Preliminary results have been presented in abstract form (Helias et al., 2007 ).

In this section we show, that the calcium peak amplitude in a spine in good approximation depends exponentially on the temporal difference of the presynaptic and the postsynaptic spiking and that the calcium influx is largest, if the presynaptic cell fires shortly before the postsynaptic cell. This makes the calcium signal an appropriate candidate carrier of information on causal correlation.

Figure 1 A illustrates the situation at a synapse subject to a spike pairing protocol. The principal mechanism causing the dependence of the postsynaptic calcium signal can readily be understood. With each presynaptic spike a small amount of glutamate is released, increasing the glutamate concentration ρ_glu at time t_glu. Glutamate binds to a fraction of the postsynaptic NMDA receptors. There is experimental evidence (Mainen et al., 1999 ) that only a fraction n(t) of the NMDA receptors bind to glutamate. n(t) reaches its maximum n₀ after a rise time of τ_rise,nmda 5 − 10 ms. We choose this maximum obtained for a single presynaptic release to be the unit of n. After the glutamate concentration ρ_glu(t) has decayed back to its resting value, the receptors unbind and n(t) decays back to 0. As long as the postsynaptic spine is not depolarized, the NMDA receptors are blocked by magnesium and thus have a low conductance (Jahr and Stevens, 1990a ,b ). Arrival of a postsynaptic back-propagating action potential (called bpAP in the following) depolarizes the spine (indicated as V_spine in Figure 1 A), removes this block and Ca²⁺ can flow into the spine. The amount of NMDA receptors opening at time t_BP therefore equals

n(t_BP) = n(t_BP − t_glu).

This is illustrated in Figure 1 B. Thus the NMDA receptor mediated conductivity depends on the relative timing Δt_syn = t_BP − t_glu. This model assumes instantaneous removal of the Mg²⁺ block. However, the detailed mechanism is more complicated. Removal of the Mg²⁺ block can happen at several time scales; a very fast one of 100 μs (Jahr and Stevens, 1990b ) and several recently found longer time scales of a few 100 ms (Kampa et al., 2004 ). The longer time scales were found to effectively narrow the time window of substantial NMDA conductance (Kampa et al., 2004 ). Here, we do not intend to capture the NMDA receptor kinetics in full detail, but rather to construct a functional yet quantitative model where the experimental results constrain the model parameters. Calcium imaging studies on spines (Nevian and Sakmann, 2004 ) show, that the timing dependence of the Ca²⁺ peak amplitude can be fitted for positive Δt_syn to a function with a rise time of τ_rise,nmda followed by a exponential decay with a single time constant of τ_nmda. These studies typically measure the timing Δt = t_post − t_pre between the spiking in the presynaptic neuron’s soma t_pre and the spiking in the postsynaptic neuron’s soma t_post (see Morrison et al., 2008 , for the definition of delays in models of synaptic plasticity). With d_BP being the delay for a bpAP and d_glu being the delay between presynaptic spike and glutamate release, the timing difference at the synapse is

In this work, we neglect the continuous binding process and instead assume that the number n of glutamate bound NMDA receptors jumps to a positive value after the rise time τ_rise,nmda, like

where H(t) denotes the Heaviside function. Assuming the bpAP to be a point event, the total Ca²⁺ influx into the spine

depends exponentially on the timing Δt where i_Ca = q₀δ(t − t_BP) is the calcium current through a single open NMDA channel during a bpAP. Our expression for q defines the effective rise time τ_rise and suggests to measure q in units of q₀n₀, i.e. we set q₀n₀ = 1. Nevian and Sakmann (2004) measured an effective rise time τ_rise = 10 ms and an exponential decay with a time constant τ_nmda = 32 ms.

The postsynaptic depolarization due to the AMPA receptor activation only causes a small NMDA conductivity, as shown in measurements of the NMDA conductivity for voltage patterns caused by spike pairing experiments (Kampa et al., 2004 ) and also directly by observing that the calcium transient in spines in spike pairing experiments is not decreased significantly by blocking AMPA receptors (Nevian and Sakmann, 2004 ). In our model, we neglect the influence of AMPA mediated depolarizations on the calcium signal.

Assuming the bpAP to be a point event is obviously an approximation. As well as the assumption, that the glutamate binding state jumps from 0 to 1 at t = t_glu + τ_rise,nmda instead of showing a continuous increase on a time scale of 5–10 ms (as found experimentally by Kampa et al., 2004 ). However, for timings Δt_syn ≥ τ_rise,nmda, it is easy to see that this does not qualitatively change the total Ca²⁺ influx (assuming i_Ca(t < 0) = 0, because

As before, q depends exponentially on the relative timing Δt_syn. The exact shape of the bpAP only enters the timing independent factor q′₀. For Δt_syn < τ_rise,nmda, the integral depends both on the shape of the resulting i_Ca(t), which in turn depends on the actual shape of the bpAP, and on the rising flank of the NMDA channel activation. A typical trace of i_Ca(t) has a half duration of 5 ms and the rise time of NMDA activation is of the same order of magnitude (5–10 ms). Hence, the slope of q(Δt_syn) in this range is large compared to the following exponential decay with τ_nmda, as confirmed experimentally (Nevian and Sakmann, 2004 ). The step like approximation (Eq. 2) therefore seems adequate. Our choice to describe the bpAP to be a point event has the consequence that in this spike pairing experiment there is no calcium influx for larger negative timings Δt_syn − τ_rise,nmda < 0. Taking into account a finite short decay time for the postsynaptic depolarization as in Shouval et al. (2002) would lead to a small calcium influx also for the post-before presynaptic timing.

The Ca²⁺ influx q leads to a transient signal which in good approximation decays in an exponential fashion (Nevian and Sakmann, 2004 , 2006 ; Waters et al., 2003 ) with a decay time of 20–200 ms (reviewed in Cavazzini et al., 2005 ). We therefore assume the calcium peak amplitude to depend linearly on the amount q of calcium influx.

So far we studied the model for the case of a spike pairing protocol, where a presynaptic action potential is paired with a postsynaptic action potential. In this section we investigate how much information about the correlation between presynaptic and postsynaptic events is contained in the Ca²⁺ signal if the spiking activity is irregular. To this end, we simulate the presynaptic and the postsynaptic spiking as Poisson processes with rates ν_i and ν_o, respectively. Both processes share a fraction of spikes that appear in the spike trains with a fixed temporal distance Δt. In the following we call these events “pair events”. Figure 2 A illustrates these processes. The strength of the temporal correlation is given by the correlation coefficient ε, which is the conditional probability of a pair event (i.e. seeing the temporally correlated presynaptic partner spike at t − Δt), given a postsynaptic spike at time t.

Let A be the set of presynaptic spike times and B the set of postsynaptic spike times within a finite interval [0, T]. The effects of the presynaptic spikes add up linearly in our model. Therefore, the number of glutamate bound NMDA receptors found at time t is

and the total Ca²⁺ influx at the point in time t_post of a postsynaptic spike is accordingly

where we again set q₀n₀ to unity. For different realizations of the presynaptic and postsynaptic spiking the elements of the set a_A,B = {q_A(t_post) | t_post ∈ B} are random variables, since A and B are random sets. We therefore need to calculate the probability density function of q_A and show that it depends on the correlation ε between pre- and postsynaptic events. Initially assume ε = 0, i.e. A and B to be two independent sets of Poisson points. In this case the points in time t_post at which q_A(t) is sampled are randomly and uniformly drawn from the interval [0, T]. Consequently the elements a ∈ a_A,B occur according to the probability density function of q_A(t), where t is a randomly and uniformly drawn point in time t ∈ [0, T]. This is the amplitude distribution ρ₀(q) of a shot noise with an exponential kernel and can be calculated recursively (Gilbert and Pollak, 1960 ) as

where γ is Euler’s constant. In the case 0 < ε ≤ 1, however, we have to distinguish two cases: Given a postsynaptic spike, with probability 1 − ε this spike is not part of a pair event. The postsynaptic event is uncorrelated with respect to the presynaptic events and thus samples q(t) at a random point in [0, T] as discussed above. So the contribution of this event to the probability density ρ(a) is (1 − ε)ρ₀(a). With probability ε the postsynaptic spike is part of a pair event. In this case, we know that there is a presynaptic spike at time t_post − Δt of which we still see the deterministic contribution to a. Since we assume the presynaptic spikes to have Poisson statistics, the existence of the presynaptic spike at t_post − Δt does not influence the statistics of the other presynaptic spikes. The latter produce a shot noise background obeying the distribution ρ₀(a′): the probability density to observe the value a = a′ + Δa equals the probability density ρ₀(a′). Hence the contribution in this case is ερ₀(a − Δa). Considering both types of postsynaptic spikes we arrive at

Figure 2 B illustrates this result. The postsynaptic events belonging to spike pairs (black bars) sample the shot noise signal at high values and cause the peak in the histogram around 0.73. Its height scales linearly with the correlation coefficient ε. The uncorrelated gray events cause a background manifested in the histogram by the large peak at low values. Its amplitude scales proportional to 1 − ε.

In conclusion, the information about correlation between the presynaptic and the postsynaptic events enters the probability distribution ρ(a) and is therefore exhibited in the amplitude distribution of Ca²⁺ transients in the postsynaptic spine: correlated postsynaptic events which come in close temporal proximity to a presynaptic spike produce a high Ca²⁺ influx.

In this section we devise a model of a postsynaptically realized counter for correlated events. The model is based on CaMKII, the most probable downstream signaling protein. CaMKII is an example of a bistable effector protein, whose transitions between its active and its inactive state are triggered by distinct calcium concentrations. Since the number of CaMKII molecules in a spine is low a statistical description is essential. Figure 3 A shows a schematic drawing of the possible transitions of a single ring molecule in our model and their dependence on the calcium concentration: If the calcium influx exceeds the highest threshold Θ_h, the number x of active CaMKII molecules is increased. We call this a plus-event. If the influx is between Θ_l and Θ_h, the amount of active molecules remains unaffected. This region is often referred to as “no man’s land” (Cormier et al., 2001 ). For a low calcium influx between Θ_b and Θ_l, the number of active molecules decreases. We call this a minus-event. If the influx is below Θ_b, the number of active molecules remains the same. Previous theoretical work (Shouval et al., 2002 ) assumes a similar dependence of the synaptic weight change on the calcium concentration, but does not take into account the intermediate concentration, where no plasticity occurs. After a high calcium event, the concentration eventually drops to levels between Θ_b and Θ_l, where it can activate the phosphatase and thus deactivates CaMKII molecules. So the increase of active molecules in our model caused by a high Ca²⁺ event is understood to be the net effect of activated minus deactivated molecules.

Cormier et al. (2001) measured the thresholds for constant calcium concentrations. It is unclear, whether these values also hold for transient calcium influx through NMDA receptors. Furthermore, for the spike pairing protocol, only relative calcium concentrations have been measured. We therefore pursue a more phenomenological approach and choose Θ_h such that a plus-event occurs, if the temporal distance between the presynaptic and the postsynaptic event Δt is in the range 0 ≤ Δt = t_post − t_pre ≤ Δt₊ = 20 ms, consistent with the potentiation window of STDP (Bi and Poo, 1998 ). Furthermore, Bi and Poo (1998) showed that the transition from LTD to LTP occurs in a relatively narrow time window symmetric around Δt = 0. Since in the present work we aim at a functional model we choose the effective rise time to be τ_rise = 0 ms, in order for the coincidence window to start at Δt = 0. Θ_h is then given by (see Figure 4 ). Subsequently, we use the absolute values of Θ_h and Θ_l of Cormier et al. (2001) to infer from the ratios and the appropriate values of Θ_l and Θ_b for our condition.

Note that this choice of thresholds has the consequence for the spike pairing protocol that for 0 ≤ Δt ≤ Δt₊ we obtain activation of CaMKII molecules, but we do not obtain deactivation for the reversed timing Δt < 0 as suggested by the experimentally observed LTD window of STDP. In the section “Sensitivity of Results to Model Assumptions” we discuss an extension of our model to incorporate the LTD window for negative relative timing and we argue that our results are invariant under this modification.

Knowing the probability distribution ρ(a) of the calcium peak amplitudes a given a postsynaptic depolarization, we can calculate the probability of occurrence p₋ = P(Θ_b ≤ a < Θ_l | postsynaptic spike) of a minus-event and p₊ = P(a ≥ Θ_h | postsynaptic spike) of a plus-event. In order to detect pairs of correlated pre- and postsynaptic events, their relative timing must satisfy Δt ≤ Δt₊. This is equivalent to the condition , meaning, that correlated events cause an influx larger than Θ_h (compare Figure 3 B).

Using Eqs. 5 and 4 for the probability of a plus-event we arrive at

Here r and C are given as in Eq. 4. An analogous calculation, which is valid under the same assumption Θ_b < Θ_l < Θ_h ≤ Δa yields

Each PSD of a spine contains a number of CaMKII molecules, which was found to be N = 80 in a typical PSD on average (Chen et al., 2005 ). We are only interested in the behavior of the two stable states of a molecule: On long time scales, a molecule can only be in the fully activated or in the completely inactive state. Let the number of active molecules be x and assume that the postsynaptic neuron spikes with a rate ν_o. Then x is a random variable and we are interested in the equilibrium distribution and its dependence on the correlation between the presynaptic and postsynaptic spike train. Each postsynaptic spike may lead to a plus-event with probability p₊(ε). In such an event, each of the N − x inactive molecules can become active. In our model, this transition happens independently for each particle with probability p. Given x active molecules, the number of molecules per time being activated is ν_o p ₊ p(N − x). Analogously, the number of molecules per time being deactivated is ν_op₋qx. Here, q is the probability of a minus-event to deactivate a particular molecules. This scenario is sketched in Figure 5 . In equilibrium, both currents must compensate, leading to the expected number of active molecules

Thus〈x〉_eq depends on the relative probabilities p₊(ε) for a plus-event and p₋(ε) for a minus-event given by Eqs. 6 and 7. Figure 5 B shows the probability distribution for the number of active molecules for an ensemble of synapses, where the presynaptic and the postsynaptic activity are correlated Poisson processes. The higher the correlation coefficient ε between presynaptic and postsynaptic activity is, the more the distribution is shifted to the right. For a derivation of the first and second moment of the probability distribution see section “Probability Distribution for the Number of Active CaMKII Molecules” in Appendix.

The amount of active CaMKII molecules is the signal that can trigger downstream processes in the postsynaptic spine. Synapse maturation caused by insertion of new AMPA receptors into the PSD is such a process. In our model we assume that the process requires the presence of a certain minimal amount X_m of active CaMKII molecules. The probability of a synapse to maturate is therefore the probability that the number x of active molecules exceeds the threshold X_m. In a premature synapse, the initial Ca²⁺ concentration is low and hence the amount of active CaMKII is low as well. We would like to know the mean time needed by the signal x, starting at x = 0, to cross the threshold X_m for the first time. This is the mean first passage time problem. We approximate the mean first passage rate by the decay rate of the slowest decaying eigenvector of the CaMKII distribution in the section “Mean First Passage Time Problem for the Number of CaMKII Molecules” in Appendix. This solution is plotted for different thresholds in Figure 6 .

We treat synapse pruning analogously. In a mature synapse the initial amount of active CaMKII is already beyond the threshold X_m. Due to pre- and postsynaptic activity, the amount x of active CaMKII may decrease or increase, depending on the rates and the pre- and postsynaptic correlation. If eventually x falls below the minimal amount X_d, the synapse dies. We choose X_d < X_m for two reasons. First, once the amount of CaMKII is high, autophosphorylation will act regeneratively (Miller et al., 2005), making the decrease of x harder. We do not model this dynamics explicitly, but rather incorporate its effect in our choice of X_d < X_m. Secondly, experiments by Le Be and Markram (2006) suggest, that there is a “period of grace” for newly formed synapses during which they are not pruned, preventing many synapses to be created and destructed in vain. Our choice is a possible implementation.

In this section we employ the correlation detection mechanism to control synaptic pruning and demonstrate its capability to regulate spike rate towards a state of low rate. Before doing so we generalize the correlation model studied in previous sections to the more realistic scenario of a temporally extended correlation between the presynaptic and the postsynaptic activity.

The correlation detection depends on the probabilities p₊ for a plus-event (high Ca²⁺ influx) and p₋ for a minus-event (low Ca²⁺ influx). A plus-event happens, whenever the influx a is in the range a ≥ Θ_h. In order to compute p₊ and analogously p₋ we need to specify the pair correlation function C_io(t₁, t₂) between the presynaptic and the postsynaptic spike train. Generally, the correlation function is defined as

Here we restrict ourselves to stationary processes, such that C_io(t₁, t₂) = C_io(t₂ − t₁) only depends on the relative timing τ = t₂ − t₁ between the presynaptic and the postsynaptic spike. Furthermore, we assume the presynaptic spikes to be Poisson events emitted at rate ν_i. The postsynaptic spikes appear with mean rate ν_o. For large τ the cross-correlation function decays to Hence we can write

C_io(τ) = ν_iν_o + c_io(τ),

where the cross-covariance c_io vanishes for large |τ|. We now know the conditional probability to observe a presynaptic spike at t_pre ∈ [t − τ, t − τ + dt] provided that a postsynaptic spike occurs at time t

For a presynaptic spike to cause a plus-event it must have appeared within the potentiation window (see Figure 4 A), i.e. 0 ≤ τ ≤ Δt₊. So the correlation detector measures the probability ε_eff exceeding chance level of finding a presynaptic spike within this potentiation window, given a postsynaptic event at t = 0. This probability is . The cross-covariance c_io(τ) decreases on a time scale comparable to the membrane time constant of the neuron, which is typically shorter than the potentiation window Δt₊ = 20 ms. An example can be seen in Figure 7 A. In this case we can make the simplification

The probabilities of plus and minus-events are then given by Eqs. 6 and 7, respectively, with ε = ε_eff . An analytic approximation for the effective correlation coefficient in the framework of linear response theory for an integrate and fire neuron model can be found in the section “Input Output Correlation of an Integrate and Fire Neuron” in Appendix. Note that the causal dependence of output spikes on the input spikes and the assumed Poisson statistics of the incoming activity leads to the input–output correlation function (Figure 7 A) which only deviates from baseline for Δt > 0. Therefore, the position of the temporal window for minus events is uncritical in this setup, as long as c_io is at baseline within this window. Thus we would obtain the same results, if the time window for minus events was at negative times.

We now have the tools to investigate synaptic pruning in the scenario depicted in Figure 7 B. A neuron initially has a number k₀ of synaptic excitatory inputs, each of which reaches the neuron via a spiny synapse with a calcium based correlation detector as described in the section “A Counter for Correlated Events” All synapses are mature and may eventually die depending on the correlation variable x: A synapse is pruned as soon as the amount x of active CaMKII undercuts a critical threshold x < X_d. The initial distribution of x over the ensemble of synapses is the eigenvector of the slowest decaying eigenmode at the initial firing frequency ν_o(0) of the neuron. The choice is justified, if we think of the initial connectivity as the outcome of a slow dynamic wiring process, during which the synaptic amount x of active CaMKII had enough time to settle in this eigenmode. In addition to the excitatory inputs, the neuron receives a static configuration of inhibitory connections. Due to the pruning process, excitatory synapses progressively die and hence the number of excitatory connections k(t) decays and the neuron’s firing rate ν_o decreases (see Figures 8 A,B). The pruning process continues until the postsynaptic neuron stops spiking. Thus, pruning is a mechanism to regulate the firing rate downwards. If the process is slow compared to the dynamics of x, we can assume the distribution of x to follow the eigenstate adiabatically. In this approximation, the development of connectivity obeys the differential equation

where η₁(ν_o, k) is the eigenvalue of the slowest decaying eigenmode. All terms are accessible: a derivation of η₁ is presented in the section “Mean First Passage Time Problem for the Number of CaMKII Molecules” in Appendix and ν_o can be calculated using Eq. 34. Therefore, Eq. 11 can be numerically integrated. Figure 8 compares this semi-analytic result with a direct simulation of the model.

The negative slope of the pruning curve in Figure 8 A is decreasing with decreasing number of incoming synapses k. The reason is the dependence of the synapse death rate on the postsynaptic firing rate: The time scale of activation and deactivation of the CaMKII molecules is determined by the postsynaptic firing rate ν_o (compare Eq. 20). Consequently, also the death rate is proportional to ν_o. This is an interesting feature, since it facilitates a firing rate homeostasis: If new synapses are created with a constant rate the input connectivity to the neuron has a stable fixed point at k* where synaptic death is just compensated by synapse creation. The firing rate of the neuron assumes a corresponding fixed firing rate ν*_o. A similar example of such a homeostasis will also be shown in the section “Synaptic Maturation and Pruning”.

In the previous section we investigated a synaptic pruning process, where all excitatory inputs are uncorrelated. However, there is evidence that coactive inputs are stabilized (Cohen-Cory, 2002 ) and therefore less likely to be pruned. Here we show, that the calcium based correlation detection mechanism naturally leads to cooperation between synapses, which stabilizes coactive inputs.

In the first setup, we investigate a neuron receiving excitatory inputs from two different sources: A pool of n_p presynaptic neurons with uncorrelated Poisson spiking activity at rate ν_i. We call these inputs the “independent inputs” in the following. The second pool of n_c Poisson spiking neurons at rate ν_i, however, generates correlated spike trains. We use the multiple interaction process (Kuhn et al., 2003 ) to produce the spike trains and call these the “correlated inputs”. The correlation coefficient 0 < c ≤ 1 is the probability of input neuron i having an input spike at time t, given neuron j has a spike at the same time. Thus, c = 0 results in uncorrelated Poisson processes, whereas for c = 1 all n_c spike trains are the same.

Initially (at t = 0) there are n_p(0) = n_c(0) = 1000 incoming synapses. Each of the synapses uses the calcium based correlation detection mechanism to determine the number of active CaMKII molecules x. Whenever x falls below the threshold X_d, the corresponding synapse dies. Figure 9 A illustrates this scenario.

Figure 9 B shows the evolution of the number of incoming connections. The synapses from independent sources exhibit a higher death rate than synapses from correlated inputs. We can readily understand this behavior: Given there is a spike at input i of the correlated pool, each of the remaining n_c − 1 inputs also delivers a spike with probability c at the same time. Thus, given a spike at the correlated input i, the expectation value for the sum of all inputs at this time is 〈w〉_mip = w(1 + c(n_c − 1)), where w is the homogeneous synaptic weight. In contrast, a spike from the uncorrelated pool only carries its own weight w < 〈w〉_mip. A higher synaptic weight results in a higher probability of the target neuron to emit a spike. Thus, the probability of the neuron to fire in response to a spike from one of the correlated inputs is higher than for an uncorrelated input. This probability is proportional to the correlation coefficient ε_eff between the presynaptic and the postsynaptic spike train (for a derivation of an analytic expression for ε_eff see Eq. 37 in the section “Correlated Poisson Input” in Appendix). With ε_eff,mip > ε_eff,Poisson, the number of active molecules x in the synapses from correlated inputs is higher than in those from independent inputs (see also Figure 5 B) and hence their death rate is lower (see also “A Counter for Correlated Events”).

Although there is no direct interaction between synapses from correlated inputs, cooperativity emerges between them and helps to stabilize these inputs in favor of the uncorrelated inputs. There is not only cooperation among synapses of the correlated pool, but also competition between the two pools: The number of uncorrelated inputs decreases with increasing correlation among the correlated inputs. This is because the synaptic deathrate increases with the firing rate of the target neuron.

The time evolution of the number of inputs n_p(t) and n_c(t) can be calculated numerically completely analogous to the section “Rate Homeostasis by Synaptic Pruning”. The system of differential equations governing the dynamics is

where η_1,Poisson(ν_o, n_p, n_c) and η_1,mip(ν_o, n_p, n_c) are the slowest decaying eigen-modes calculated according to the section “Mean First Passage time Problem for the Number of CaMKII Molecules” in Appendix.

In this section, we extend the scenario of Figure 9 A by incorporating a process of synapse creation. Again we have two pools of input neurons: The first pool of N_p neurons has Poisson spiking activity with rate ν_i, the second pool of N_c neurons has Poisson activity with pair correlation c and the same rate ν_i. Excitatory synapses from both pools can exist in either of two states, premature or mature as depicted in Figure 10 A. Newly created synapses are in the premature state, lacking AMPA receptors. Their synaptic weight is 0. A synapse becomes mature, if the number of active molecules x exceeds a threshold X_m. The mature synapse has the synaptic weight w > 0. This synapse dies, if the number of active molecules falls below a threshold X_d. Initially there are no premature synapses n_p,pre(0) = n_c,pre(0) = 0, and there are as many mature synapses from independent inputs as from correlated inputs n_p,mat(0) = n_c,mat(0) > 0. Premature synapses are constantly created: For each of the N_p − n_p,mat − n_p,pre presently unestablished connections from independent sources, the rate of realization is λ_pre. Premature synapses from correlated sources are created analogously with the same rate λ_pre.

The evolution of the connectivity in the presence of maturation and pruning is shown in Figure 10 B. The connectivity approaches an equilibrium state after t 600 s. The number of mature synapses from correlated inputs n_c,mat increases, while the number of synapses from uncorrelated inputs n_p,mat decreases (upper two traces). The two values saturate at different levels, such that n_c,mat > n_p,mat. Initially, the number of premature synapses increases for both input pools (lower two traces). In equilibrium, the numbers of premature synapses saturate at different levels n_p,pre > n_c,pre. The explanation for this observation is the same as in the previous subsection: Correlatedly activated synapses exhibit cooperation, they experience a higher correlation coefficient between the input spike train and the output spike train produced by the neuron. Hence, for synapses from the correlated pool, the maturation rate of premature synapses is higher and the deathrate of mature synapses is lower as compared to the uncorrelated inputs. As in the section “Synaptic Pruning”, the target neuron becomes selective for the correlated pool.

A similar differential equation as Eq. 12 quantitatively describes the evolution of the connectivity in this model. Its numerical solution (Figure 10 B, gray) corresponds well to a direct simulation of the system (black). The observed deviations are due to synapse maturation: Each synapse entering the mature state has a number x of active CaMKII just above threshold X_m. This perturbs the equilibrium distribution of x for the mature synapses and hence influences their death rate. The direction of this influence depends on the relative position of X_m with respect to the equilibrium 〈x〉_eq: For X_m > 〈x〉_eq the observed death rate is lower than the analytic estimate and vice versa.

In the present work we describe a novel model of the synaptic mechanisms controlling synapse pruning and synapse maturation. To our knowledge, this is the first model of structural plasticity based on the microscopic dynamics of the single synapse. Instead of constructing a phenomenological model, we use recent experimental findings to identify plausible postsynaptic mechanisms and to constrain the choice of model parameters. We analyze the cortical scenario of irregular spiking activity and quantify the dependence of structural plasticity on the correlation between the presynaptic and the postsynaptic activity. Functionally, the dynamics acts as a Hebbian learning rule for synaptogenesis and pruning. Our model can be understood as a biological plausible implementation of a cascade synapse (Rubin and Fusi, 2007 ) for spiking activity, where the number x of active CaMKII molecules plays the role of internal states. However, in contrast to the cascade model, once x falls below the critical threshold the synapse is lost and cannot be reactivated by subsequent potentiating events. Here we provide a full statistical treatment of the number of active molecules. This is essential, since due to the low number of molecules (N 80) fluctuations cannot be neglected. Previous work on STDP (Cai et al., 2007 ; Shouval and Kalantzis, 2005 ) demonstrated, that stochastic fluctuations can change the phenomenology of a learning rule.

We show that there are two distinct rate regimes: In the low rate regime up to ν_i 10 Hz, our model works as a correlation detector in the sense that synchronously activated synapses on the same dendrite stabilize. This frequency limit is a direct consequence of the NMDA unbinding time constant τ_nmda. Hence the model constitutes a viable mechanism for a network using temporal (correlation) codes at low rates. At higher rates (ν_i > 10 Hz), synapse stabilization occurs irrespective of correlations between input and output. This is in accordance with recent experimental evidence that higher activity leads to increased synapse formation (Le Be and Markram, 2006 ) and to classical LTP induced by tetanic stimulation protocols. Furthermore, we show that by controlling synaptic pruning our model exhibits several desirable features for a neuronal network: The synaptic pruning rate increases with the postsynaptic firing rate. This enables the homeostasis of firing rate while synapses are continuously created, a finding obtained earlier using abstract rate based models of structural plasticity (van Ooyen et al., 1995 ). Moreover, the synapses targeting the same neuron naturally exhibit cooperation and competition. These emerging phenomena render the proposed microscopic mechanism relevant for the theory of learning in neuronal networks: The evolution of connectivity is sensitive to correlations in the inputs and hence neurons become selective for coactive inputs. This stabilization of coactive inputs has been proposed earlier based on experimental evidence, reviewed in (Cohen-Cory, 2002 ). Our treatment explains, how cooperation and competition are mediated solely by the identical postsynaptic activity experienced by synapses. Previous theoretical work (Kempter et al., 1999 , 2001 ) showed cooperativity to emerge from the interplay of spike timing dependent learning rules with the spiking dynamics in the framework of spike response models. In general, a quantitative understanding of the interaction between a spike timing based plasticity rule and the integrate-and-fire dynamics is a hard problem. Here we provide such an analysis for our specific learning rule assuming all-to-all spike interaction. The analysis allows us to obtain semi-analytic expressions for the evolution of network structure. Specifically, we present results for the case of incoming irregular Poisson activity as well as for the case of correlated inputs generated by a multiple-interaction-process (Kuhn et al., 2003 ).

For analytical convenience we use an integrate and fire neuron model with δ-shaped postsynaptic currents, as commonly used in network simulations and theoretical works. These currents can cause an immediate spiking response to an incoming spike. More realistic neuron models have postsynaptic currents with finite rise times and hence also the input–output correlation function shows a finite latency of τ_response 5 ms in response to an input spike. Since we aim at a consistent theory of the interaction between learning rule and neuronal dynamics, we compensate for the lack of latency by reducing the measured glutamate binding time τ_rise 5…10 ms to τ_rise = 0 ms. As in the natural setting the whole mass of the input–output correlation function falls into the time window τ₊ of the learning rule. A mismatch of the time constants would not change the observed phenomena qualitatively, but slightly reduce the sensitivity to correlations. We can extend our analysis to more realistic and more complicated neuron models, if an expression for the corresponding integral input–output correlation is known.

The probabilities p and q for activation and deactivation of CaMKII molecules by large and small calcium events respectively, are set to p = q = 0.01 without experimental reference. Once experimental data are available, these parameters need to be reconsidered. The smaller the values, the narrower the distribution of active CaMKII molecules. This increases the sensitivity to correlation but also the latency of the distribution in following transient changes in the correlation. The threshold X_d for the minimum number of active CaMKII molecules required for survival can be chosen such that the synapse turnover in a neuronal network reaches experimentally observed values of 7% of all synapses per week (Stettler et al., 2006 ). As soon as an experimental value for X_d is available we can check whether our model consistently relates the two experimental measures. We use a multiple interaction process to generate incoming spiking activity with higher order correlations. The choice is motivated by analytical convenience rather than by biological realism. However, we expect other models of higher order correlation to exhibit qualitatively similar results for cooperation and competition between synapses.

Due to our assumption of a point event like bpAP (see also “Spike Time Dependence of Postsynaptic Calcium Concentration”) there is no calcium influx through NMDA receptors, if the postsynaptic depolarization precedes the presynaptic release of glutamate. For our model the consequence is that post-before-pre pairings do not cause deactivation of CaMKII molecules (see also “A Counter for Correlated Events”). Hence if the intention was to explain STDP, we would not be able to reproduce the part of the rule expressing LTD. Experimentally, Nevian and Sakmann (2004) observe calcium influx through NMDA receptors only for the pre-before-post condition. Additionally, there is influx through voltage dependent calcium channels (VDCC) opened by the bpAP. Also AMPA receptors depolarize the spine and hence lead to calcium entry through coactivated NMDA receptors and VDCCs. However, the resulting NMDA conductivity is small (Kampa et al., 2004 ) and calcium transients in spines are not significantly decreased by blocking AMPA receptors (Nevian and Sakmann, 2004 ). So here we assume the main pathway for calcium influx into the spine to be NMDA receptors and we neglect the additional sources of calcium. As long as their contribution to the calcium signal is small, they constitute a background elevation of the overall calcium level in the spine and we obtain qualitatively similar results. If, however, their contributions were of comparable size as the calcium transient through NMDA receptors, a temporally close post-before-pre timing could lead to sufficient calcium influx to deactivate CaMKII and hence create a LTD window. Previous biophysical models of STDP (Shouval et al., 2002 ) assumed a finite decay time for the postsynaptic depolarization. This as well leads to an intermediate calcium influx for the post-before-pre condition. Our model can be extended in the same way. However, as argued in the section “Rate Homeostasis by Synaptic Pruning”, for the results presented here, the exact temporal position of the LTD window is irrelevant. The observed rate homeostasis only depends on the proportionality between the synapse death rate and the postsynaptic firing rate, resulting from the CaMKII dynamics. The cooperativity between synapses (see “Cooperation and Competition by Spatial Input Correlation”) only requires a Hebbian type plasticity, i.e. plus events are caused by a postsynaptic spike following a presynaptic spike in close succession.

The treatment of the CaMKII dynamics as a Markov process implicitly assumes the high (plus) and low (minus) calcium events to occur uncorrelatedly with the respective probabilities p₊ or p₋. However, the occurrence of a plus event or a minus event depends on the current value n(t_post), the fraction of glutamate bound NMDA receptors at the time of the postsynaptic spike. n(t) has a time constant τ_nmda = 32 ms. Hence, if postsynaptic spikes occur with arbitrarily small inter-spike-intervals, the probability to observe a plus event is slightly higher after a previous plus event. Nevertheless, for realistic postsynaptic spike trains, small inter-spike-intervals are rare due to refractoriness, so the neglect of the temporal correlation of n(t) is well justified (see also section “Mean First Passage Time Problem for the Number of CaMKII Molecules” in Appendix, Figure 12 B). A more thorough treatment must take into account the actual auto-correlation function of the postsynaptic spiking activity. The dynamics of the discrete amount of active CaMKII molecules is mapped to a continuous system. By comparing the continuous analytic probability density to the discrete numerical solution of the probability mass function, we found that this approximation is sufficiently accurate as long as the total number of molecules is large enough (N > 30). Furthermore, we approximate the CaMKII dynamics as a diffusion with a x-independent diffusion constant, whereas the original problem leads to a x-dependent diffusion term. Again, the nearly perfect agreement of the numerical solution (taking into account x-dependent diffusion) with the analytically obtained Gaussian, proves this approximation adequate (see Figure 11 ). However, if the distribution approaches the saturation limits x = 0 or x = N, we observe pronounced deviations. This may occur for highly correlated spiking and at excessive firing rates. The results for the synaptic death and maturation rates are based on the slowest decaying eigenvalue η₁ of the CaMKII activation distribution. This component has the largest time constant. If the ensemble of synapses is initialized with a distribution containing contributions different from the slowest decaying eigenvector, we observe transient deviations in the pruning rate. These transients decay typically 15 times faster than the slowest eigenvector So after a sufficiently long time, the CaMKII distribution obeys the analytical solution in good approximation. For calculating the equilibrium state of connectivity, the slowest decaying eigenvector is the exact description if the number x of active molecules of a new synapse is drawn from this eigenvector. But even for non-stationary connectivity, the distribution of activated CaMKII across synapses follows the structure quasi adiabatically and the approximation is well fulfilled as shown by comparing analytical results with direct simulations (see “Rate Homeostasis by Synaptic Pruning” and “Cooperation and Competition by Spatial Input Correlation”). The larger deviations in Figure 10 are attributed to this approximation: Each synapse entering the mature state has a number x of active CaMKII just above threshold X_m. This perturbs the equilibrium distribution, explaining the deviation of the simulated from the analytically obtained connectivity structure.

Recent experiments (Stettler et al., 2006 ) demonstrated synaptogenesis and synaptic death to occur even in the adult cortex. However, their functional relevance for the neural network still remains to be illuminated. Theoretical considerations (Chklovskii et al., 2004 ; Stepanyants et al., 2002 ) suggest rewiring in networks to provide the dominating substrate of information storage. In artificial neural networks optimization of wiring was already shown to contribute significantly to memory capacity (Knoblauch, 2006 ; Knoblauch et al., 2007 ). More generally, restructuring of connectivity allows biological systems to optimize their circuitry to fulfill a specific set of functions. Understanding the mechanisms of wiring optimization will also be beneficial for technical systems (e.g. integrated circuits), since due to their essentially two dimensional nature, the number of contacts existing at any point in time is a precious and limited resource. For the above purposes a phenomenological model of structural plasticity is sufficient. However, a microscopic model is required to uncover the underlying biological mechanisms and resulting limitations. Moreover, an understanding of the control mechanisms of connectivity may contribute to the development of medical protocols to promote plasticity after neural lesions, e.g. as experienced after a stroke. How a system can exhibit a plastic structure and yet acquire and maintain its functionality is still a matter of research. To answer this question it is not sufficient to investigate synaptogenesis and synaptic death under conditions of stationary activity, but we rather have to consider the close interplay between the structural dynamics and the correlation dynamics exposed by the present work. The framework presented here provides the tools for this endeavor; we are now in the position to investigate structural plasticity in recurrent neural networks.

Suppose there is a finite amount of CaMKII in the PSD. Each of the N molecules can either be active or inactive. Transitions between these states are triggered by calcium influx, where a high calcium amplitude a ≥ Θ_h causes each inactive molecule to be activated with probability p, and a low calcium amplitude Θ_b ≤ a < Θ_l causes each active molecule to be inactivated with probability q. The rate of high calcium events is λ₊ = ν_o p₊(ε), the rate of low calcium events λ₋ = ν_o p₋(ε), where ν_o denotes the postsynaptic firing rate. Given x active molecules, the number of molecules being activated per time is λ₊p(N − x), the current of molecules being deactivated is λ₋qx. Since these rates depend on the state x of the system, the process has the Markov property. Thus, the system is uniquely defined by the transition probability P(x, y) from state x into state y. Let λ = ν_o(p₊(ε) + p₋(ε)) be the rate of events that change the state of the system, the probability that the event was a plus-event and be defined analogously. Given an event (either plus or minus), the transition probability is

where is the binomial distribution and H the Heaviside function. We are interested in the equilibrium distribution ρ(x), which must fulfill

To obtain the non-trivial solution, we have to take into account the normalization condition This results in the inhomogeneous linear system of N + 1 equations

which can be solved numerically. Since the transition probability P(x, y) is a positive stochastic matrix, i.e. it fulfills P(x, y) ≥ 0 and according to the Perron-Frobenius theorem (MacCluer, 2000 ) its largest eigenvalue is 1 and the respective eigenvector ρ(x) is unique with positive entries. This guarantees a unique solution of Eq. 14 with the desired properties of a probability distribution.

In order to obtain information about the equilibrium density ρ(x) we investigate its first and second moments. Suppose, x obeys the distribution ρ(x, t) at time t with a well defined first and second moment. In this case we can determine the infinitesimal time evolution of μ(t) = 〈x(t) 〉 and σ²(t) =〈(x(t) −〈x〉)²〉 for short times Δt. For the mean we obtain

where the identity (mean of the binomial distribution) was used. The first moment of the distribution fulfills the differential equation

In equilibrium, the mean value is

In addition, we can calculate its temporal evolution from Eq. 15 to be

So the time scale on which the distribution approaches its equilibrium is given by

Similarly, the infinitesimal time evolution of the variance unfolds to

where we used and k² = Np(1 − p) + (Np)². We obtain the differential equation

In equilibrium, we determine

Figure 11 compares the numerical solution of Eq. 14 to a Gaussian with mean μ_eq and standard deviation σ_eq, demonstrating that the Gaussian approximation is sufficient.

In the model the number of active CaMKII molecules x determines whether a synapse maturates or dies. Both decisions are triggered by the crossing of different thresholds. Thus we have to calculate the mean time until x passes the threshold X for the first time. This is known as the mean first passage time problem with an absorbing boundary at x = X. Here we are interested in an approximation for the mean first passage time, or equivalently, for the threshold passing rate.

Figure 5 shows that the equilibrium distribution ρ(x) in absence of any absorbing boundary is well described by a Gaussian distribution. This observation suggests a mapping of the discrete first passage time problem to a continuous one with and the corresponding replacement of the discrete probability distribution ρ by a continuous probability density function

In the following we omit the tilde; x and ρ name the continuous variables. The time evolution of the probability density ρ(x) can be described by the Fokker–Planck equation

where S(x) denotes the probability flux operator. An explicit form can be determined under the following assumptions: (1) The process can be described as a diffusion. (2) The equilibrium distribution of the process is given by a Gaussian where μ_eq and σ_eq are determined by Eqs. 16 and 17. The stationary probability density has to fulfill so from Eq. 19 it follows that S(x)ρ₀(x) = S₀ = const. The constant must be S₀ = 0, since the probability current vanishes at it vanishes for all x. We intend to describe the process as a diffusion. Therefore, S should contain only first derivatives in x and the flux operator must have ρ₀ as its stationary solution, leading to

The constant D controls the time scale of the process. We determine D such that the dynamics of the diffusion process Eq. 19 matches the dynamics of our process for the first moment of ρ. In doing so, we follow Ricciardi et al. (1999) and determine the infinitesimal drift term A₁(x) from the master equation of the process, which yields

Subsequently we require A₁(x) = −D(x − μ_eq), which fixes

To simplify the notation, we rescale the variable The flux operator now reads A separation ansatz for the probability density ρ(y, t) = φ(y)e^ηt turns the Fokker–Planck equation 19 into the corresponding eigenvalue problem

To find the eigenvalues of L_FP, we follow Risken (1996) and transform Eq. 21 into a Hermitian operator, using the following transformation

which satisfies

Hence the transformed operator L is

Let φ(y) be an eigenvector of L_FP, then obviously ψ = U⁻¹(y)φ(y) is an eigenvector of L with the same eigenvalue, so we have to solve

This is the differential equation of the harmonic oscillator known from quantum mechanics. The general solution of this differential equation is a linear combination of confluent hypergeometric functions (Abramowitz and Stegun, 1974 )

The last equality follows from M(a, b, z) = e^zM(b − a, b, −z) where the constant c must be determined using the normalization condition The boundary conditions constrain the choice of possible values for η as well as for the constants r. Without loss of generality, we assume that the threshold X has to be crossed from above (synapse pruning). The case in which the threshold represents an upper boundary (synapse maturation) can be handled analogously by interchanging the roles of plus and minus-events and the associated quantities. At the upper end of the interval, the probability flux has to vanish, i.e. this is a totally reflecting boundary

The lower end of the interval is given by the threshold X. This is an absorbing boundary, but does not require the probability density to vanish, since the rate of minus-events is limited and hence the threshold passing rate remains finite even if φ(X) > 0. The probability for the system in state x to cross the threshold X after having received a minus-event is

where the last expression is the cumulative density function of the binomial distribution. Figure 12 shows the typical steep increase of the exit probability in the vicinity of the threshold X.

The flux η through the boundary can be calculated as

The probability density can be approximated by a linear Taylor series here, since p_exit(x) vanishes except in a small range near the threshold (see Figure 12 ). Using the flux operator S, the flux through the boundary can be expressed as

Taken together we obtain the boundary condition

Using the explicit form of ρ given by Eq. 22, and its derivative

where we exploit the identity the boundary condition Eq. 23 assumes the form

which fixes the constant r as

Finally, we simultaneously solve Eqs. 28 and 27 for η by a simple numerical bisection method. By construction, the solutions η₀ > η₁>… are eigenvalues of the Fokker–Planck operator Eq. 21. Since eventually the system has to pass the threshold, there cannot be an equilibrium eigenstate with η₀ = 0; all eigenstates are decaying. We are interested in the largest negative eigenvalue, which determines the slowest decaying eigenvector Eq. 22, i.e. the largest η₁ < 0 which solves Eq. 21 and fulfills the boundary conditions. Being initially (t = 0) in this eigenstate, the probability of finding the system still in the interval [0, X] after time t is

where η₁ is called the escape rate or the inverse mean first passage time. The eigenvector in Eq. 22 resulting from the considerations above is plotted in Figure 6 together with results of direct simulations. The decay rate expected from Eq. 29 agrees well with the simulation. In this analytical treatment we have assumed that plus and minus-events occur independently of the history of previous events with probabilities and respectively. However, this is an approximation, since plus and minus-events depend on the signal n(t), which has the time scale τ_nmda. Thus, given a postsynaptic spike at t_post,1 and given n(t_post,1) > Θ_h, a plus-event is produced. n(t) then decays exponentially with the time constant τ_nmda or it jumps to even higher values, if a presynaptic spike arrives. Therefore a later postsynaptic spike t_post,2 > t_post,1 definitely generates another plus-event as long as i.e. within a finite time window after the first event. Consequently, the probability for a plus-event to follow another plus-event is slightly higher than The same is true for minus-events. But since a neuron’s spike train exhibits refractoriness and the probability is low that the second spike occurs within the time window determined by τ_nmda, the correlated occurrence of plus (or analogously minus) events has only a small impact. This is verified in Figure 12 B, where the probability of k consecutive plus-events is displayed. For events caused by the spike train of an integrate and fire neuron, the curve is near its theoretically expected value while for spikes with Poisson statistics, the correlated occurrence of plus-events can be explained by with

Here we calculate the correlation coefficient of the incoming spike train at a given excitatory synapse and the output spike train of an integrate and fire neuron, where the correlation coefficient is defined by Eq. 10. If we have incoming Poisson spikes at a stationary rate ν_i, the neuron fires in a stationary fashion with rate ν_o as well. We can then rewrite the cross-correlation function as C_io(Δt) = ν_o(Δt | given an input spike at t = 0)ν_i and formally the correlation coefficient in Eq. 10 reads

For the particular case of Δ-shaped postsynaptic currents we find an explicit expression for Eq. 30 as follows. The membrane potential V is governed by

with reset potential V_r and spike threshold V_th. Here t_i is the time of the i-th incoming event and w_i its weight. For each event, the weight w_i is a random variable with distribution K. This allows us to calculate Eq. 30 also for the case where a neuron receives N correlated spike trains (see section “Correlated Poisson Input”). The incoming events are assumed to be Poisson events with rate ν. Then the membrane potential is a Markov process with distribution P(V, t) which obeys the evolution equation

We closely follow the derivation of the Fokker–Planck equation for the integrate and fire neuron given by Gerstner and Kistler (2002) , but generalize it for arbitrary weight distributions K. As we assume stationary input, we set t′ = 0 without loss of generality. For short times Δt the transition probability becomes

and the membrane potential distribution after time Δt is

Taking the limit Δt → 0 yields

A Kramers–Moyal expansion can be made under the assumption that the weight w is small and the rate ν is high. This is also known as the diffusion limit:

In the limit of small 〈w〉_K only the first and the second moment of K enter the Fokker–Planck equation. This determines the drift constant μ and the diffusion constant σ² to

In equilibrium, the firing rate ν_o of the neuron can be calculated by the formula (Brunel, 2000 ; Brunel and Hakim, 1999 ; Ricciardi, 1977 )

Finally, we employ linear response theory to calculate the input–output correlation coefficient in Eq. 30 for one particular incoming synapse in the limit of small weights w. First assume, that every event has the same weight w. Since the incoming spike rate is stationary, without loss of generality, we can assume an incoming spike to occur at time t = 0. This is a small perturbation of the membrane potential and shifts the distribution by weight w to the right, which can be taken into account in the Fokker–Planck equation 32 as

To first order in w, the incoming spike can be described as a Δ shaped perturbation of the mean input μ(t) = μ + τwΔ(t). Since the events are Poisson, their auto-correlation is flat, thus in the instant t = 0⁻ before the spike arrives the neuron’s membrane potential obeys the equilibrium distribution. The effect of the perturbation on the firing rate can be treated in the framework of linear response theory for small w << V_th − μ. The response of the firing rate can then reads

where ν_o is the equilibrium firing rate and h(t) is the impulse response of the firing rate with respect to a Δ-shaped perturbation of the mean input. So Eq. 30 becomes

where is the step-response of the firing rate for t → ∞ with respect to a unit-step change of the mean input. Up to linear order in w

In the second step we employ Eq. 34 and in the third summarize terms by the quantity ω(μ, σ) called “dc-susceptibility” in linear response theory. Using Eq. 36 in Eq. 35 results in the desired explicit expression. A related approach can be found in De la Rocha et al. (2007) .

The input–output correlation coefficient can also be calculated in the case, where the neuron receives spikes from N excitatory inputs of weight w generated by a multiple-interaction-process (mip) Kuhn et al. (2003) with mother spike rate ν_mother and copy probability c. In addition, the neuron receives uncorrelated excitatory and inhibitory Poisson input which we treat as a Gaussian white noise background input with mean μ_bg and standard deviation σ_bg. As in “Poission Input” we are interested in the correlation coefficient in Eq. 30 of a given input i ∈ [1,…,N] and the output spike train. A spike at time t_i at input i, indicates that there is a mother-spike that was successfully copied to input i. Hence, with probability c each of the remaining N − 1 inputs also has a spike at time t_i. Therefore the expectation value of the weight of the composed event is 〈w〉_mip = (1 + (N − 1)c)w. For small 〈w〉_mip <<V_th − μ the we can apply the framework of linear response theory as before and calculate Eq. 30 using Eqs. 36 and 35 as

where

The last two equations using μ_mip and result from Eq. 33.

A synapse connecting an axon of neuron i to the dendrite of neuron j is stored in a list associated with neuron i. Without loss of generality we restrict the discussion to the case of a single synapse between i and j. In a distributed simulation, each process stores only the part of the list referring to the resident target neurons (see Morrison et al., 2005 ). The parameters and dynamic variables stored for each synapse j are: the synaptic weight w_j, the synaptic delay composed of the dendritic delay and the axonal delay the number x_j of active CaMKII molecules, and the boolean variable mature_j indicating whether the synapse is mature (True) or premature (False). Axonal delays and the NMDA receptor rise time can be taken into account as long as (compare Morrison et al., 2005). However, in our simulations, we assume for sake of simplicity. We also assume the time constant τ_nmda to be identical for all synapses. In this case the fraction of glutamate bound NMDA receptors n(t) as defined in Eq. 3 is the same for all axonal synapses of neuron i. In addition, each neuron stores the time t_old of its last spike. The spike distribution algorithm invokes the method spike(t) (see Algorithm 1) for each spike t of the presynaptic neuron i. It propagates the dynamics of each synapse in the local target list. Here we use the function T(x, n′) (see Algorithm 2) as an abbreviation for the transitions of x due to the number n′ = n(t_post) of glutamate bound NMDA receptors as defined in the section “A Counter for Correlated Events”. To evolve the synaptic dynamics, we need the spike history of the postsynaptic neuron. Here we closely follow Morrison et al. (2007) and use the same definition of update_register(t) and get_history(t₁, t₂). An extension, however, is necessary, because in the presence of structural plasticity the number of incoming connections of a neuron may change over time: The access counter counter_SP of each entry in the spike register has to be adapted whenever a synapse is created or dies. We define the method start_recording(t_first) (see Algorithm 3) called on the postsynaptic neuron when a new synapse is created. The argument is the earliest the new synapse will start accessing the spike history (excluding t_first). Analogously, stop_recording(t_last) is called when a synapse dies (see Algorithm 4). Its argument is the latest point in the history (including t_last) considered by the synapse before dawn (see Algorithm 1).

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

Many thanks to Birgit Kriener for very helpful discussions and to Tom Tetzlaff for comments on an earlier version of the manuscript. Ursula Körner and Andreas Knoblauch contributed essential inspiring ideas in the initial phase of the project. The implementation of the synapse model is the result of constructive discussions with the NEST initiative (www.nest-initiative.org). Partially funded by DIP F1.2, BMBF Grant 01GQ0420 to the Bernstein Center for Computational Neuroscience Freiburg, and EU Grant 15879 (FACETS), and “The Next-Generation Integrated Simulation of Living Matter” project, part of the Development and Use of the Next-Generation Supercomputer Project of the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.