{'fraction_non_alphanumeric': 0.029450384543237666, 'fraction_numerical': 0.01026341881716108, 'mean_word_length': 4.5747684493576335, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structural Plasticity Controlled by Calcium Based Correlation Detection

Hebbian learning in cortical networks during development and adulthood relies on the presence of a mechanism to detect correlation between the presynaptic and the postsynaptic spiking activity. Recently, the calcium concentration in spines was experimentally shown to be a correlation sensitive signal with the necessary properties: it is confined to the spine volume, it depends on the relative timing of pre- and postsynaptic action potentials, and it is independent of the spine's location along the dendrite. NMDA receptors are a candidate mediator for the correlation dependent calcium signal. Here, we present a quantitative model of correlation detection in synapses based on the calcium influx through NMDA receptors under realistic conditions of irregular pre- and postsynaptic spiking activity with pairwise correlation. Our analytical framework captures the interaction of the learning rule and the correlation dynamics of the neurons. We find that a simple thresholding mechanism can act as a sensitive and reliable correlation detector at physiological firing rates. Furthermore, the mechanism is sensitive to correlation among afferent synapses by cooperation and competition. In our model this mechanism controls synapse formation and elimination. We explain how synapse elimination leads to firing rate homeostasis and show that the connectivity structure is shaped by the correlations between neighboring inputs.


INTRODUCTION
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 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: 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 infl ux 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 confi rmed 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 fl uctuations 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 affi nity to calcium than CaMKII. Therefore, they become active at lower calcium concentrations and counteract the phosphorylation of CaMKII. This is in agreement with the fi nding 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 simplifi ed neuron models (Butz et al., 2008;Dammasch et al., 1986) or plasticity rules depending on the fi ring 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 fl uctuations 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 infl uence 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 fi ring 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 fi ndings 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 fi ndings. In the section "Ca 2+ 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  using the computationally effi cient 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).

SPIKE TIME DEPENDENCE OF POSTSYNAPTIC CALCIUM CONCENTRATION
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 infl ux is largest, if the presynaptic cell fi res shortly before the postsynaptic cell. This makes the calcium signal an appropriate candidate carrier of information on causal correlation. Figure 1A 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 0 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 1A), removes this block and Ca 2+ can fl ow 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 1B. Thus the NMDA receptor mediated conductivity depends on the relative timing Δt syn = t BP − t glu . FIGURE 1 | (A) Sketch of an experiment measuring the calcium concentration in spines during a spike pairing protocol. The calcium concentration inside the spine is measured while varying the temporal difference Δt syn between the onset t glu of the glutamate concentration and the arrival of a postsynaptic bpAP t BP at the spine. (B) Rise in glutamate concentration at t glu in the synaptic cleft causes NMDA receptors in the spine to bind to glutamate. The fraction n (t) of glutamate bound NMDARs jumps to its maximum (n 0 = 1 for simplicity) after time τ rise,nmda . As long as the spine is not depolarized, the NMDA pore is blocked by Mg 2+ . A postsynaptic spikes causes a bpAP to arrive at the spine at t BP and to unblock the currently glutamate bound NMDA receptors n(t BP ), so Ca 2+ can enter the spine. The total amount of Ca 2+ infl ux q is proportional to n(t BP ) and therefore depends on the relative timing Δt syn . This model assumes instantaneous removal of the Mg 2+ block. However, the detailed mechanism is more complicated. Removal of the Mg 2+ 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 2+ peak amplitude can be fi tted 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 defi nition 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 2+ infl ux into the spine (2) depends exponentially on the timing Δt where i Ca = q 0 δ(t − t BP ) is the calcium current through a single open NMDA channel during a bpAP. Our expression for q defi nes the effective rise time τ rise and suggests to measure q in units of q 0 n 0 , i.e. we set q 0 n 0 = 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 signifi cantly by blocking AMPA receptors (Nevian and Sakmann, 2004). In our model, we neglect the infl uence 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 2+ infl ux (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′ 0 . 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 fl ank 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 confi rmed 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 infl ux for larger negative timings Δt syn − τ rise,nmda < 0. Taking into account a fi nite short decay time for the postsynaptic depolarization as in Shouval et al. (2002) would lead to a small calcium infl ux also for the postbefore presynaptic timing.
The Ca 2+ infl ux q leads to a transient signal which in good approximation decays in an exponential fashion 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 infl ux.

Ca 2+ TRANSIENTS CAUSED BY CORRELATED IRREGULAR SPIKING
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 2+ 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 fi xed temporal distance Δt. In the following we call these events "pair events". Figure 2A illustrates these processes. The strength of the temporal correlation is given by the correlation coeffi cient ε, 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. where we again set q 0 n 0 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 ρ 0 (q) of a shot noise with an exponential kernel and can be calculated recursively (Gilbert and Pollak, 1960) as , γ = .
0 577215665 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 − ε)ρ 0 (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 Δ = − Δ − a t e rise nmda ( ) / τ τ to a. Since we assume the presynaptic spikes to have Poisson statistics, the existence of the presynaptic spike at t post − Δt does not infl uence the statistics of the other presynaptic spikes. The latter produce a shot noise background obeying the distribution ρ 0 (a′): the probability density to observe the value a = a′ + Δa equals the probability density ρ 0 (a′). Hence the contribution in this case is ερ 0 (a − Δa). Considering both types of postsynaptic spikes we arrive at with e rise nmda Figure 2B 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 coeffi cient ε. 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 2+ transients in the postsynaptic spine: correlated postsynaptic events which come in close temporal proximity to a presynaptic spike produce a high Ca 2+ infl ux. Ӎ 0 73 (with τ rise = 5 ms). The black dots result from a realization of the stochastic process (T = 10 5 s) using Scientifi c Python (Jones et al., 2001), the gray curve illustrates (Eq. 5). The peak amplitude is proportional to the correlation coeffi cient ε = 0.1. Uncorrelated postsynaptic events cause the peak at low values of a.

A COUNTER FOR CORRELATED EVENTS
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 3A 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 infl ux exceeds the highest threshold Θ h , the number x of active CaMKII molecules is increased. We call this a plus-event. If the infl ux 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 infl ux between Θ b and Θ l , the number of active molecules decreases. We call this a minus-event. If the infl ux 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 2+ 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 infl ux 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 . 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 modifi cation.
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 Δ := −Δ a t h e nmda / τ ≥ Θ , meaning, that correlated events cause an infl ux larger than Θ h (compare Figure 3B).
Using Eqs. 5 and 4 for the probability of a plus-event we arrive at x increases x decreases ρ(a) a FIGURE 3 | (A) CaMKII molecules exhibit bistability: An inactive state and a highly phosphorylated, active state. Transitions between the two states are triggered by Ca 2+ infl ux. If the infl ux exceeds a threshold Θ h , the molecule switches from inactive to active. Low Ca 2+ infl ux between Θ b and Θ l deactivates the molecule. (B) At each postsynaptic spike, Ca 2+ can enter the spine and can change the amount x of active CaMKII molecules. The probability density function ρ(a) of the Ca 2+ infl ux at the points in time of a postsynaptic spike a = q(t post ) determines the direction of the effect on x (same data as in Figure 2B). According to this effect, it can be divided into different regions: For Θ b ≤ a < Θ l; the number of active molecules is reduced. This event occurs with probability p − (Eq. 7). For a ≥ Θ h the number is increased, occurring which probability p + (Eq. 6). In the regions between Θ l and Θ h and below Θ b there is no infl uence on x.

Frontiers in Computational
Here r and C are given as in Eq. 4. An analogous calculation, which is valid under the same assumption 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 ν o p − 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 plusevent and p − (ε) for a minus-event given by Eqs. 6 and 7. Figure 5B 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 coeffi cient ε between presynaptic and postsynaptic activity is, the more the distribution is shifted to the right. For a derivation of the fi rst 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 2+ 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 fi rst time. This is the mean fi rst passage time problem. We approximate the mean fi rst 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.

RATE HOMEOSTASIS BY SYNAPTIC PRUNING
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 2+ infl ux) and p − for a minus-event (low Ca 2+ infl ux). A plus-event happens, whenever the infl ux 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 1 , t 2 ) between the presynaptic and the postsynaptic spike train. Generally, the correlation function is defi ned as .
Here we restrict ourselves to stationary processes, such that C io (t 1 , t 2 ) = C io (t 2 − t 1 ) only depends on the relative timing τ = t 2 − t 1 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 lim ( ) .
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 4A), i.e. 0 ≤ τ ≤ Δt + . So the correlation detector measures the probability 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 7A. In this case we can make the simplifi cation 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 coeffi cient in the framework of linear response theory for an integrate and fi re 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 7A) 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 7B. A neuron initially has a number k 0 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 fi ring frequency ν o (0) of the neuron. The choice is justifi ed, 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 confi guration 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 fi ring rate ν o decreases (see Figures 8A,B). The pruning process continues until the postsynaptic neuron stops spiking. Thus, pruning is a mechanism to regulate the fi ring 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.

FIGURE 7 | (A)
Input-output correlation function normalized to the postsynaptic spike rate. The peak drops to baseline (ν i = 5 Hz, gray) on a time scale which is shorter than the width of the potentiation window Δt + = 20 ms. The area below the peak is the effective correlation coeffi cient ε eff (see Eq. 10). (B) A neuron N receiving k Poisson spiking inputs of rate ν i = 5 Hz via synapses s 1 …s k . Each synapse measures the correlation between its input spike train and the spiking activity of the neuron N. A synapse is pruned as soon as its number of active CaMKII has fallen below a threshold x < X d = 30. In this approximation, the development of connectivity obeys the differential equation where η 1 (ν o , k) is the eigenvalue of the slowest decaying eigenmode. All terms are accessible: a derivation of η 1 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 8A is decreasing with decreasing number of incoming synapses k. The reason is the dependence of the synapse death rate on the postsynaptic fi ring rate: The time scale of activation and deactivation of the CaMKII molecules is determined by the postsynaptic fi ring rate ν o (compare Eq. 20). Consequently, also the death rate is proportional to ν o . This is an interesting feature, since it facilitates a fi ring rate homeostasis: If new synapses are created with a constant rate the input connectivity to the neuron has a stable fi xed point at k* where synaptic death is just compensated by synapse creation. The fi ring rate of the neuron assumes a corresponding fi xed fi ring rate ν* o . A similar example of such a homeostasis will also be shown in the section "Synaptic Maturation and Pruning".

COOPERATION AND COMPETITION BY SPATIAL INPUT CORRELATION
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.

SYNAPTIC PRUNING
In the fi rst 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 coeffi cient 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 9A illustrates this scenario. Figure 9B 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 fi re 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 coeffi cient ε 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 5B) 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 fi ring 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.

SYNAPTIC MATURATION AND PRUNING
In this section, we extend the scenario of Figure 9A by incorporating a process of synapse creation. Again we have two pools of input neurons: The fi rst 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 10A. 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 10B. 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 coeffi cient 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 10B, 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 infl uences their death rate. The direction of this infl uence 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.

DISCUSSION
In the present work we describe a novel model of the synaptic mechanisms controlling synapse pruning and synapse maturation. 0 05 1 s synapses per second in otherwise the same input scenario as in Figure 9A. Initially, there are n p,mat (0) = n c,mat (0) = 3540 mature excitatory synapses from each input pool and n p,pre (0) = n c,pre (0) = 0 premature synapses. Additionally, there are inputs from 1120 inhibitory synapses not subject to pruning. The traces show (top to bottom): mature synapses from correlated pool n c,mat , mature synapses from uncorrelated pool n p,mat , premature synapses from uncorrelated pool n p,pre , premature synapses from correlated pool n c,pre . Black: simulation, gray: analytical approximation.

Helias et al. Calcium based correlation detection
To our knowledge, this is the fi rst model of structural plasticity based on the microscopic dynamics of the single synapse. Instead of constructing a phenomenological model, we use recent experimental fi ndings 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) fl uctuations cannot be neglected. Previous work on STDP (Cai et al., 2007;Shouval and Kalantzis, 2005) demonstrated, that stochastic fl uctuations 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 fi ring rate. This enables the homeostasis of fi ring rate while synapses are continuously created, a fi nding 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(Kempter et al., , 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-andfi re dynamics is a hard problem. Here we provide such an analysis for our specifi c learning rule assuming all-to-all spike interaction. The analysis allows us to obtain semi-analytic expressions for the evolution of network structure. Specifi cally, 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).

SENSITIVITY OF RESULTS TO MODEL ASSUMPTIONS
For analytical convenience we use an integrate and fi re 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 fi nite rise times and hence also the input-output correlation function shows a fi nite 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 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 infl ux 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 infl ux through NMDA receptors only for the prebefore-post condition. Additionally, there is infl ux 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 signifi cantly decreased by blocking AMPA receptors (Nevian and Sakmann, 2004). So here we assume the main pathway for calcium infl ux 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 suffi cient calcium infl ux to deactivate CaMKII and hence create a LTD window. Previous biophysical models of STDP (Shouval et al., 2002) assumed a fi nite decay time for the postsynaptic depolarization. This as well leads to an intermediate calcium infl ux 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 fi ring 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.

SCOPE AND ACCURACY OF THE ANALYTICAL TREATMENT
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 justifi ed (see also section "Mean First Passage Time Problem for the Number of CaMKII Molecules" in Appendix, Figure 12B). 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 suffi ciently 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 fi ring rates.
The results for the synaptic death and maturation rates are based on the slowest decaying eigenvalue η 1 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 suffi ciently 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 fulfi lled 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.

FUNCTIONAL ROLE OF STRUCTURAL PLASTICITY
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 artifi cial neural networks optimization of wiring was already shown to contribute signifi cantly to memory capacity (Knoblauch, 2006;Knoblauch et al., 2007). More generally, restructuring of connectivity allows biological systems to optimize their circuitry to fulfi ll a specifi c set of functions. Understanding the mechanisms of wiring optimization will also be benefi cial 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 suffi cient. 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 suffi cient 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.

Numeric solution of the equilibrium distribution
Suppose there is a fi nite amount of CaMKII in the PSD. Each of the N molecules can either be active or inactive. Transitions between Calcium based correlation detection these states are triggered by calcium infl ux, 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 fi ring 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 defi ned 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, + which can be solved numerically. Since the transition probability P(x, y) is a positive stochastic matrix, i.e. it fulfi lls P(x, y) ≥ 0 and ∀ : , = = x P x y y N Σ 0 1 ( ) , 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.

Analytic approximation of the equilibrium distribution
In order to obtain information about the equilibrium density ρ(x) we investigate its fi rst and second moments. Suppose, x obeys the distribution ρ(x, t) at time t with a well defi ned fi rst and second moment. In this case we can determine the infi nitesimal time evolution of μ(t) = 〈x(t)〉 and σ 2 (t) = 〈(x(t) − 〈x〉) 2 〉 for short times Δt. ( ) (mean of the binomial distribution) was used. The fi rst moment of the distribution fulfi lls 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 infi nitesimal time evolution of the variance unfolds to σ ρ λ eq eq e where μ eq and σ eq are determined by Eqs. 16 and 17. The stationary probability density has to fulfi ll ∂ρ ∂t = 0, so from Eq. 19 it follows that S(x)ρ 0 (x) = S 0 = const. The constant must be S 0 = 0, since the probability current vanishes at x = − 1 2 it vanishes for all x. We intend to describe the process as a diffusion. Therefore, S should contain only fi rst derivatives in x and the fl ux operator must have ρ 0 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 fi rst moment of ρ. In doing so, we follow Ricciardi et al. (1999) and determine the infi nitesimal drift term A 1 (x) from the master equation of the process, which yields (17) 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 suffi cient.

Mean fi rst passage time problem for the number of CaMKII molecules
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 fi rst time. This is known as the mean fi rst passage time problem with an absorbing boundary at x = X. Here we are interested in an approximation for the mean fi rst 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 fi rst 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 December 2008 | Volume 2 | Article 7 | 16

Helias et al. Calcium based correlation detection
To fi nd the eigenvalues of L FP , we follow Risken (1996) and transform Eq. 21 into a Hermitian operator, using the following transformation ( ) 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 confl uent hypergeometric functions (Abramowitz and Stegun, 1974) The last equality follows from M(a, b, z) = e z M(b − a, b, −z) where the constant c must be determined using the normalization condition 1 = ∫ φdx. 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 minusevents and the associated quantities. At the upper end x N = + 1 2 of the interval, the probability fl ux has to vanish, i.e. this is a totally refl ecting 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 fi nite 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 fl ux η 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 fl ux operator S, the fl ux through the boundary can be expressed as Taken together we obtain the boundary condition  where η 1 is called the escape rate or the inverse mean fi rst 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 + p and − p , 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 defi nitely generates another plusevent as long as n t i For the particular case of δ-shaped postsynaptic currents we fi nd 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 fi re 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 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 fi rst and the second moment of K enter the Fokker-Planck equation. This determines the drift constant μ and the diffusion constant σ 2 to μ = τν〈w〉 K σ 2 = τν〈w 2 〉 K .
In equilibrium, the fi ring rate ν o of the neuron can be calculated by the formula (Brunel, 2000;Brunel and Hakim, 1999;Ricciardi, 1977) 1 Finally, we employ linear response theory to calculate the input-output correlation coeffi cient 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 fi rst 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 fl at, 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 fi ring rate can be treated in the framework of linear response theory for small w << V th − μ. The response of the fi ring rate can then reads ν o (t | input spike at t = 0) = ν o + wτh(t) t ≥ 0, where ν o is the equilibrium fi ring rate and h(t) is the impulse response of the fi ring rate with respect to a δ-shaped perturbation of the mean input. So Eq. 30 becomes where ∫ ′ ′ ∞ 0 h t t ( )d is the step-response of the fi ring 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).

Correlated poisson input
The input-output correlation coeffi cient 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 coeffi cient 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 The last two equations using μ mip and σ mip 2 result from Eq. 33.

ALGORITHMIC IMPLEMENTATION OF SYNAPSE MATURATION AND SYNAPSE DEATH
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 d d d  Morrison et al., 2005). However, in our simulations, we assume d j A + = τ rise,nmda 0 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 defi ned 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 defi ned in the section "A Counter for Correlated Events". To evolve the synaptic dynamics, we need the spike history of the postsynaptic neuron. Here we Algorithm 1 spike(t) Require: t old last presynaptic spike processed, n = n(t old ) Ensure: t old = t, n = n(t) on exit for all postsynaptic neurons j do if mature j then send spike (  x ← x − Binomial(x, q) end if return x Algorithm 3 start_recording(t fi rst ) start at beginning of spike register while t SP <= t fi rst do counter SP ← counter SP + 1 move to next element end while