text:
Article | Open | Published:

# Associative properties of structural plasticity based on firing rate homeostasis in recurrent neuronal networks

## Abstract

Correlation-based Hebbian plasticity is thought to shape neuronal connectivity during development and learning, whereas homeostatic plasticity would stabilize network activity. Here we investigate another, new aspect of this dichotomy: Can Hebbian associative properties also emerge as a network effect from a plasticity rule based on homeostatic principles on the neuronal level? To address this question, we simulated a recurrent network of leaky integrate-and-fire neurons, in which excitatory connections are subject to a structural plasticity rule based on firing rate homeostasis. We show that a subgroup of neurons develop stronger within-group connectivity as a consequence of receiving stronger external stimulation. In an experimentally well-documented scenario we show that feature specific connectivity, similar to what has been observed in rodent visual cortex, can emerge from such a plasticity rule. The experience-dependent structural changes triggered by stimulation are long-lasting and decay only slowly when the neurons are exposed again to unspecific external inputs.

## Introduction

Network plasticity involves connectivity changes at different levels. Changes in the strength of already existing synapses are known as functional plasticity, whereas structural changes of axonal or dendritic morphology, as well as the creation of new and deletion of already existing synapses, is known as structural plasticity. During certain stages of development, axons and dendrites have been shown to grow and degenerate depending on neuronal activation1,2,3. Structural changes, however, are not limited to the extent and shape of neurites, but also include more subtle alterations in spines and boutons. Spine remodeling on excitatory cells due to neuronal activity has been observed in vitro using organotypic hippocampal cultures4,5,6 and in vivo in sensory and in motor cortex of rodents7,8,9,10,11. Not only changes on postsynaptic spines, but also changes in presynaptic structures have been reported in organotypic hippocampal slice cultures12,13,14.

The exact rules governing activity-dependent structural changes, however, are still not understood. Computational models have tried to shed light on this issue by simulations, showing for example that many observed features of cortical connectivity could be achieved through the interaction of multiple plasticity mechanisms15,16, and that homeostatic regulation of neuronal activity on multiple time scales is necessary in order to stabilize Hebbian changes17,18. Homeostatic plasticity, in this context, usually refers to a regulation of neuronal connectivity that result in stabilization of neuronal activity at a set point. There is growing evidence for homeostatic regulation of cortical connectivity (for a review, see Turrigiano19). More recently, homeostatic regulation of cortical activity has been demonstrated in vivo in rodent visual cortex20,21,22,23. Hebbian plasticity, on the other hand, is used to describe mechanisms that change connections between two neurons based on the correlation between their respective activities. Many of these aspects were discussed at a recent conference devoted to the interaction of Hebbian and homeostatic plasticity24.

But is it the case that the associative principles defining Hebbian learning must rely on pre-post correlations that are available exclusively at the level of individual synapses? Or could associative learning in networks also emerge from plasticity rules that are based on homeostatic principles on the level of whole neurons? The latter has been first proposed by Dammasch25 and, to our knowledge, has not been followed up since then. He proposed that Hebbian learning could emerge as a network property by using an algorithm based on firing rate homeostasis of individual neurons, with no reference to correlation. The idea that associative learning could also emerge from the principle of homeostasis brings an important new aspect into the discussion of integrating Hebbian and homeostatic plasticity. Experimental data usually describe the effects of plasticity on connectivity, but there are still many details missing. Therefore, it is currently not an easy task to differentiate between three scenarios how associative learning might arise: It could emerge from a correlation-based Hebbian learning rule, it might arise as network effect of a learning rule based on homeostasis, or it could occur as a combination of both.

In this paper, we will explore these questions by testing the idea proposed by Dammasch25, using a reimplementation of his algorithm in a more modern modeling framework. We use a structural plasticity rule based on firing rate homeostasis recently implemented in NEST26,27 and a recurrent network of leaky integrate-and-fire (LIF) excitatory and inhibitory neurons. We show that a strongly interconnected assembly of neurons emerges when these neurons are jointly stimulated with stronger external input. We then test the associative properties in an experimentally well-documented scenario, employing a simple model for the maturation of circuits in the primary visual cortex of rodents (V1).

It has been shown that neurons in adult V1 have an increased probability to be synaptically linked to other neurons that have a similar preference for visual features28. Later, it was demonstrated that this feature-specific bias in connectivity was not present at the time of eye opening, and it developed only after some weeks of visual experience29, suggesting that plastic mechanisms would shape the maturation of V1 circuits through visual experience. Surprisingly, a feature-specific bias in connectivity was also shown to develop after eye opening in dark-reared mice lacking any visual input30. The idea of activity dependent plastic mechanisms shaping the maturation of these networks, however, was not ruled out. Spontaneous retinal activity is also present on dark-reared animals31, and plasticity could refine V1 networks based on patterned activity received by pairs of neurons that share some of these inputs. Moreover, the relationship between connectivity and neuronal response to natural movies was not as strong for dark-reared mice as for normally reared mice30, suggesting that plastic mechanisms contribute to the maturation of V1 circuits in an essential way, and the full maturation of feature specific connectivity would depend on visual experience anyway.

Recently, Sadeh et al.32 have shown that a bias for feature specific connectivity can emerge in balanced random networks of LIF neurons with a synaptic plasticity rule that combines Hebbian and homeostatic mechanisms. They used, however, a functional plasticity rule, which limits a direct comparison to experimental data on connectivity. Here, we show that feature-specific connectivity can also emerge in a network of LIF neurons with a structural plasticity rule in which correlations are implicitly evaluated through the random combination of presynaptic and postsynaptic elements in the network, and which does not require that synapses keep track of the presynaptic activity. Moreover, we observed long-lasting structural after-effects of stimulation. This property is compatible with the notion of a persistent memory, which is not in every moment reflected by activity.

## Methods

### Network simulations

We simulate a recurrent network of N = 12500 current-based LIF neurons, of which N E  = 0.8 N are excitatory and N I  = 0.2 N are inhibitory. The sub-threshold dynamics of the membrane potential V i of neuron i obeys the differential equation

$${\tau }_{m}\frac{d{V}_{i}}{dt}=-{V}_{i}+{\tau }_{m}\sum _{j}\,{J}_{ij}\,{s}_{j}(t-d),$$
(1)

where τ m is the membrane time constant. The synaptic weight J ij from a presynaptic neuron j to a postsynaptic neuron i is the peak amplitude of the postsynaptic potential and depends on the type of the presynaptic neuron. Excitatory connections have a strength of J E  = J = 0.1 mV. Inhibitory connections are stronger by a factor g = 8 such that J I  = −gJ = −0.8 mV. A spike train $${s}_{j}(t)={\sum }_{k}\delta (t-{t}_{j}^{k})$$ consists of all spikes produced by neuron j. Some of these neurons represent the external input. The lumped spike train of all external neurons to a given neuron in the network is modeled as a Poisson process of rate ν ext , and external input to different neurons is assumed to be independent. All synapses have a constant transmission delay of d = 1.5 ms. When the membrane potential reaches the firing threshold Vth = 20 mV, the neuron emits a spike to all postsynaptic neurons and its membrane potential is reset to V r  = 10 mV and held there for a refractory period of tref = 2 ms.

The indegree is fixed at 0.1N I for inhibitory to inhibitory and inhibitory to excitatory connections, and at 0.1N E for excitatory to inhibitory synapses. After connections of these types are established they remain unchanged throughout the simulation. In contrast, excitatory to excitatory (EE) connections are initially absent and emerge only from the structural plasticity rule. All simulations were conducted using the NEST simulator33,34]. Numerical values of all parameters of our model are again collected in Table 1.

### Homeostatic structural plasticity (SP)

The SP model used in our work has been recently implemented in NEST27. The implementation combines precursor models by Dammasch35, van Ooyen & van Pelt36 and van Ooyen37. This model has been employed before to study the rewiring of networks after lesion or stroke26,38,39, the specific properties of small-world networks40, the emergence of critical dynamics in developing neuronal networks41, and neurogenesis in the adult dentate gyrus42,43. All these models, however, included a distance-dependent kernel for the formation of new synapses, which is not part of the NEST implementation27 that was used in our present study.

#### Neuronal activity and synaptic elements

The EE connections in the network are volatile and undergo continuous remodeling, controlled by the SP algorithm. In its first versions, the model had continuous representations of pre- and postsynaptic densities, which were used for accessing connectivity between two neurons42. Later, it was adapted to have a discrete number of axonal and dendritic elements, which are combined to form synapses between neurons38. In the original model, the electrical activity of a neuron is represented by its intracellular calcium concentration, which is a lowpass filtered version of its time-dependent firing rate. In this paper, we use the lowpass filtered spike train of neuron i as a measure for its instantaneous firing rate r i

$${\tau }_{r}\frac{d{r}_{i}(t)}{dt}+{r}_{i}={s}_{i}(t\mathrm{).}$$
(2)

The time constant of the lowpass filter was chosen as τ r  = 10 s throughout all our simulations.

Excitatory neurons are assigned a target rate ρ and a set of pre- and postsynaptic elements, which can be interpreted as axonal boutons and dendritic spines, respectively. The number of presynaptic elements, zpre, and postsynaptic elements, zpost, evolves in dependence of the neuron’s firing rate according to

$$\frac{d{z}_{i}^{k}(t)}{dt}=-\frac{1}{{\beta }^{k}}({r}_{i}(t)-\rho ),k\in \{{\rm{pre}},\,{\rm{post}}\},$$
(3)

where i is the index of the neuron, and βk is a growth parameter. We use a target rate ρ = 8 Hz and growth parameter β = 2 for both pre- and postsynaptic elements of all excitatory neurons, unless stated otherwise.

Previous work on SP has used different functions to describe how these elements change with the neuron’s activity, such as linear38,41,42,43, gaussian26,27,39, or logistic40. Since there is currently no direct experimental data showing how the number of these elements vary with firing rate, we chose a generic linear function to implement a simple phenomenological model of firing rate homeostasis. We use the same growth rule and the same parameters for both types of elements of all neurons in our network. In the original model26, free elements that are not engaged in a synapse will decay with a certain rate. In the model considered here, however, free elements do not decay with time. We did run test simulations considering the decay, and found that our main results were not altered (see Supplementary material).

At regular intervals Δt = 100 ms, the structural plasticity rules are applied to delete already existing and create new synaptic contacts. Numerical values of the parameters are summarized in Table 2.

#### Synapse deletion and creation

At regular intervals Δt, when rewiring is scheduled, a neuron may have more or less pre- and post-synaptic elements than it has actual synapses. In that case, synapses are either deleted or created, in order to match the number of elements to the number of active synaptic contacts. At fixed intervals Δt, the number of postsynaptic (presynaptic) elements is compared to the number of existing incoming (outgoing) synapses of each neuron

$${\rm{\Delta }}{z}_{j}^{{\rm{p}}{\rm{r}}{\rm{e}}}={z}_{j}^{{\rm{p}}{\rm{r}}{\rm{e}}}-\sum _{i}{C}_{ij}\quad {\rm{a}}{\rm{n}}{\rm{d}}\quad {\rm{\Delta }}{z}_{i}^{{\rm{p}}{\rm{o}}{\rm{s}}{\rm{t}}}={z}_{i}^{{\rm{p}}{\rm{o}}{\rm{s}}{\rm{t}}}-\sum _{j}{C}_{ij},$$
(4)

where C is the matrix containing the number of synapses between presynaptic neurons j and postsynaptic neurons i. If the neuron has more synaptic contacts than synaptic elements (Δzk < 0), synapses are deleted. The $$|{\rm{\Delta }}{z}^{k}|$$ synapses to be deleted are randomly chosen among the existing contacts a neuron has. After a synapse has been deleted due to a loss of presynaptic (postsynaptic) elements, the corresponding postsynaptic (presynaptic) elements remain available for a new connection. If the neuron has more elements than contacts ($${\rm{\Delta }}{z}^{k}\mathrm{ > }0$$), the neuron is considered to have Δzk free synaptic elements. All free synaptic elements in the network are randomly combined into pairs of pre- and postsynaptic elements to form new synapses. The number of synapses formed is limited by both the total number of pre- and the total number of postsynaptic elements in the full network. Each newly created synapse has a fixed strength J. Multiple synapses between the same pair of neurons are allowed, but auto-synapses (self-connection of a neuron onto itself) are not. See Butz & van Ooyen26 and Diaz-Pier et al.27 for more details on the implementation of the model.

### Subgroup stimulation

The networks were first grown without structured input, with all excitatory neurons receiving external Poisson input with the same rate νext. Stimulation was started after 750 s, when enough EE connections were grown, and, apart from small fluctuations, all neurons fired at their target rate. During the stimulation period, a subgroup comprising 10% of the excitatory neurons received an increased external input (1.1νext for 150 s). After stimulation, the external input was set back to its original value (νext) for all excitatory neurons. Both the activity and the connectivity of the network were monitored for 5500 s.

### Visual stimulation protocol

For the visual cortex simulations, we consider a network similar to the one described in Sect. Network Simulations. The stimulation only starts after the networks have created enough EE connections such that the preset firing rate can be maintained. Visual stimulation is simulated by providing the excitatory neurons with Poisson input the rate of which depends on the orientation of the stimulus, see Sadeh et al.32 for details of the protocol. The baseline firing rate νext is the same as during the growth period, the modulation depends on the orientation of the visual stimulation θ and the preferred orientation (PO) θPO of the input and a modulation gain parameter μ

$${\nu }_{{\rm{mod}}}(\theta ,{\theta }_{{\rm{PO}}})={\nu }_{{\rm{ext}}}[1+\mu \,\cos \,\mathrm{(2(}\theta -{\theta }_{{\rm{PO}}}))]\mathrm{.}$$
(5)

Each neuron is assigned a parameter θPO, which is randomly drawn from a uniform distribution on [0°, 180°). During the stimulation phase, a different θ is randomly drawn from a uniform distribution on the interval [0°, 180°) and presented to all excitatory neurons for a duration of tst = 1 s. We use a modulation μ = 0.15, and the stimulation protocol consists of presenting a total of Nst = 5000 different stimuli. After the stimulation, the external input to all excitatory neurons is once again set to its initial non-modulated value of νext, and the network is simulated for another tpost = 10000 s. Numerical values of all parameters regarding the stimulation protocol are collected in Table 3.

### Spike train analysis

The spike count correlation between a pair of neurons i and j was calculated as the Pearson correlation coefficient

$${R}_{ij}=\frac{{c}_{ij}}{\sqrt{{c}_{ii}\,{c}_{jj}}},$$
(6)

where c ij is the covariance between spike counts extracted from spike trains x i and x j of two neurons, and c ii is the variance of spike counts extracted from x i . Correlations were calculated from spike trains comprising 20 s of activity, using bins of size 10 ms.

The irregularity of spike train of neuron i was measured as the coefficient of variation of its inter-spike intervals

$${{\rm{CV}}}_{i}=\frac{{\sigma }_{i}}{{\mu }_{i}},$$
(7)

where μ i is the mean and σ i is the standard deviation of the inter-spike intervals extracted from the spike train of neuron i (duration 20 s).

### Data availability

The datasets generated during and analysed during the current study are available from the corresponding author on reasonable request.

## Results

We start by growing recurrent networks of excitatory and inhibitory LIF neurons. All synaptic connections are static, except for EE connections, which are initially absent and grow according to a structural plasticity rule that implements firing rate homeostasis (see Methods for more details). Before stimulation, we characterized the networks formed under the influence of the structural plasticity rule for uniform (untuned) external stimulation.

### Grown networks are random in absence of structured input

The target rate is set to ρ = 8 Hz for all excitatory neurons, which is the expected firing rate of excitatory neurons for the parameter set we are using and 10% EE connectivity. As expected, average EE in- and outdegree increase from 0 until stabilizing at approximately 1000, corresponding to an average EE connectivity of 10% (Fig. 1B and C). The plasticity rule is always active and individual EE connections are still being created and deleted, but we consider the network to be in equilibrium at this point, when excitatory neurons fire on average at their target rate.

Since multiple synapses between the same pair of neurons are allowed, we also looked into the distribution of the number of synapses between pairs of pre- and post-synaptic neurons. Figure 1D shows that this distribution is roughly a Poisson distribution. If individual contacts can be considered as independent random variables with a Poisson distribution, the in- and outdegree of individual neurons would also follow a Poisson distribution and satisfy μ ≈ σ2. Figure 1C shows, however, that the distribution of in- and outdegrees have σ2 < μ. In the SP model, in- and outdegree distributions change for different distribution of target rates. We also performed simulations in which target rates were drawn from broader distributions, yielding also broader distributions of in- and outdegree,and the main results of feature specific connectivity were not altered (see Supplementary material). Since a thorough study of the effect of target rate on degree distributions was beyond the scope of this paper, for simplicity we fixed ρ = 8 Hz for all excitatory neurons.

For the parameters used on our simulations, an inhibition dominated random network of LIF neurons with 10% connection probability has been shown to have low firing rates, as well as asynchronous and irregular (AI) spike trains44. Figure 1E shows the activity of the network in equilibrium, after 750 s of simulation. The population raster plot (Fig. 1E) indicates the network is in an AI state. Pairs of neurons fire with low correlation coefficient (CC, bin size 10 ms) (Fig. 1F), and individual spike trains have a coefficient of variation (CV) around 0.7 (Fig. 1F). As expected from the choice of ρ, the firing rate of excitatory neurons has a mean value of 8 Hz in equilibrium (Fig. 1F). Firing rates, pairwise correlation and irregularity were calculated from spike trains recorded during 20 s of simulation.

The networks grown according to the SP rule do exhibit some non-random features (see Supplementary Material). However, we consider these deviations from random networks as small. After 750 s of simulation, networks with a target rate ρ = 8 Hz for all excitatory neurons, and a growth parameter β = 2 for all synaptic elements have an essentially random structure and an activity that can be classified as AI.

### Time constant of growth process depends on target rate and the growth parameter for synaptic elements

We then determined the time scale of network growth. To that end we simulated a network in which all neurons receive Poisson input with the same rate νext until the equilibrium was reached. We considered the network to be in equilibrium when connectivity (i.e. in- and outdegree distributions) are stable. In this state, individual synapses are still plastic and are created, deleted and recreated as the simulation runs. For these simulations, all neurons had the same target firing rate ρ, but we simulated networks with different values set for ρ.

The evolution of EE connectivity is plotted in Fig. 1B. The time scale of that growth process depends on how fast synaptic elements grow (parameter β), but also on the target rate ρ (β = 2 and ρ = 8 Hz on Fig. 1B). The number of new synapses created in the network depends on the number of available free elements. This number, in turn, depends on the growth parameter β for the synaptic elements. Therefore, it comes as no surprise that the network growth also depends on β. The observed relation between the time constant β of network growth and target rate ρ, however, is not self-explaining, as for networks with fixed incoming inhibitory connections and fixed external input, the gain in firing rate depends non-linearly on the actual number of incoming excitatory connections to excitatory neurons44.

To better understand this dependence, we defined a growth time constant τgrowth = s/α, where s is the plateau value of the connectivity, calculated here as the average connectivity for the last 4 discrete time points in a long enough simulation, and α is the highest (mostly initial) slope extracted from the connectivity time series. We then plotted τgrowth against the growth parameter β of synaptic elements, for different values of ρ (Fig. 2A). The range of parameters considered here spanned more than 3 orders of magnitude, so we used a log-log plot to represent it. The exponent, however, is 1, indicating a linear dependency between β and τgrowth. Also, the growth of EE connectivity is faster for simulations with higher ρ, as can be seen in Fig. 2A. The time constant of EE connectivity growth is, therefore, a function of both β and ρ. In order to express this dependency of the growth time constant on ρ,we extracted the coefficient γ = τgrowth/β from the simulated data (Fig. 2B). Finally we rescaled the time axis of Fig. 1B by dividing it by τgrowth for given values of β and ρ, obtaining a growth curve in time units of γ β (Fig. 2C).

We found that the growth process is stable throughout at least 3 orders of magnitude of the parameter β. An even slower process would lead to exceedingly long simulation times, but we have no reason to believe that this would destabilize the system. A value of β could easily be chosen such that network growth would happen in hours or days, matching experimental data of structural plasticity. Making the process faster, however, may destabilize the system, as the rate of new contacts created would increase. Apart from mathematical constraints, there are also biological limits regarding speed, of course, as a fast system also requires high turnover rates and efficient transport of proteins and other molecules.

From now on, all simulated networks have a target rate ρ = 8 Hz for all excitatory neurons, all synaptic elements have a growth parameter β = 2.

### Associative properties of SP

In his paper, Dammasch25 suggested that the compensation algorithm, on which this particular SP model26,27 is based, could implement a form of Hebbian learning. This becomes clear when analyzing the equations by Butz & van Ooyen26 describing the expected change in connectivity induced by the SP algorithm:

$${\rm{\Delta }}{C}_{i,j}=(\begin{array}{ll}{C}_{i,j}\frac{{\rm{\Delta }}{z}_{j}^{pre}}{{z}_{j}^{pre}} & {\rm{for}}\,{\rm{\Delta }}{z}_{j}^{pre}\,\mathrm{ < }\,0\\ {C}_{i,j}\frac{{\rm{\Delta }}{z}_{i}^{post}}{{z}_{i}^{post}} & {\rm{for}}\,{\rm{\Delta }}{z}_{i}^{post}\,\mathrm{ < }\,0\\ \frac{{\rm{\Delta }}{z}_{j}^{pre}{\rm{\Delta }}{z}_{i}^{post}}{max(\sum _{k}{\rm{\Delta }}{z}_{k}^{pre},\sum _{k}{\rm{\Delta }}{z}_{k}^{post})} & {\rm{for}}\,{\rm{\Delta }}{z}_{j}^{pre}\,\mathrm{ > }\,0\,{\rm{and}}\,{\rm{\Delta }}{z}_{i}^{post}\,\mathrm{ > }\,0\end{array}\mathrm{.}$$
(8)

Please note that we have adapted Eq. 8 from26 to match the nomenclature used in the present paper. Also, we only account for the simplest case here, in which only EE connections are plastic, and where the formation of synapses does not depend on distance. Eq. 8 clearly shows that the connectivity increase happens when two neurons are simultaneously in a low activity state, implementing Hebbian learning through a covariance rule. On the other hand, in periods of elevated activity, the connectivity decreases in an unspecific manner, similar to a weight-dependent synaptic scaling, affecting both pre- and postsynaptic elements.

The associative properties of SP are tested by stimulating a subgroup of neurons and quantifying the changes in their connectivity. After the network was grown and all excitatory neurons were firing roughly at their target rate, a subgroup of the excitatory neurons was stimulated with higher external input for a duration of 150 s. Upon stimulation onset, the instantaneous firing rate of the stimulated subgroup increases (Fig. 3C), pushing the neurons away from their target rate and triggering a rewiring of their dendrites and axons. The connectivity of the neurons then decays (Fig. 3E) until the instantaneous rate again reaches its target value, and the average connectivity then stabilizes. When specific stimulation stops, and all excitatory neurons receive again un-tuned external input, the firing rate of the subgroup drops to a level below the set point, as the excitatory amplification by the recurrent network is now reduced due to the deletion of connections during stimulation. The firing rate below the set point triggers once again rewiring, and the neurons create new pre- and postsynaptic elements. Since the effect of specific stimulation on the firing rate of the remaining network is much smaller (Fig. 3C and E), there are more synaptic elements available from the subgroup, and it is more likely for them to create connections within the subgroup than with neurons outside of it (Fig. 3B).

The structural plasticity rule is continuously remodeling the network and there are still changes in connectivity even after the firing rate of the subgroup reaches its target (approximately 300 s after the end of specific stimulation, see Fig. 3C). In this case, however, all neurons are firing on average at their target rate, and rewiring is slower, as the change in the number of elements is proportional to r(t) − ρ. The higher connectivity that was formed within the stimulated subgroup lasts for a longer period after activity is back to what it was before stimulation (Fig. 3C and E).

### SP leads to feature-specific connectivity in a simple model for the maturation of V1

We were then interested in testing the associative properties in a biologically more realistic scenario. As an example, we consider a simple model for the maturation of V1, in which different stimuli are presented consecutively, and neurons respond to them according to their own functional preferences. Once the network was formed and is in equilibrium, excitatory neurons were driven by external input that was tuned to stimulus orientation to simulate visual experience. Each neuron received external input as a Poissonian spike train, the rate of which was modulated according to a tuning curve with an input PO (θPO) that was randomly assigned at the beginning of the simulation. More specifically, the modulation depended on the difference between θPO and the stimulus orientation (SO), which changed randomly at fixed time intervals. According to the stimulation protocol, neurons receive a slightly higher external input when the presented SO is similar to their input PO, which also entails a higher output rate (Fig. 4A–C).

Before oriented stimulation, the connectivity between excitatory neurons was random, and in particular, there was no bias of connectivity with regard to the PO of neurons (Fig. 4D and E left columns). After the presentation of 5000 different stimuli, however, the connectivity pattern of the excitatory neurons change, and neurons are more likely to connect to neurons with similar PO (Fig. 4D middles column). The connectivity is, therefore, modulated according to the difference in the POs of pairs of neurons (Fig. 4E middle column). During modulated stimulation, neurons fire either at a lower or at a higher rate than their target rate, according to the difference between the SO and their PO. Whenever they fire lower than their target rate, they create synaptic elements. In contrast, if they fire higher than their target rate, they delete synaptic elements. Neurons with similar PO increase their number of elements at the same stimulus periods and, therefore, have a higher probability of creating synapses between each other.

Figure 4D and E refer to average connectivity between classes of neurons, and synapses contribute the same to connectivity irrespective of whether they are established between two different pairs of neurons or between the same pair. Figure 4F show how the number of synapses between pairs of connected neurons are modulated according to the difference in their respective POs. After visual stimulation, neurons are more likely to create multiple synapses to other neurons with similar PO, as compared to neurons with different PO. Another aspect to notice on Fig. 4D and E is the conspicuous drop in average connectivity following the feature-specific stimulation. One possible explanation for this drop is the non-linearity of the input-output response curve of LIF neurons. During oriented stimulation the average output is higher, although the average input to the neurons is still the same as before stimulation. We performed simulations with a static network with 10% connectivity and the same modulated external input, mimicking visual experience. We found that the average firing rate of the neurons increased by 1 Hz during stimulation (see Supplementary material). In the network with SP such an increase in firing rate would lead to a decrease in average connectivity to keep average firing rate at its target value.

### Structural changes induced by stimulation decay very slowly

As in the case of subgroup stimulation (Fig. 3E), the network structure formed during visual stimulation persists long after the end of stimulation. After 5000 stimuli, we set the external input to all excitatory neurons back to a uniform value of νext to see what happens to the SP induced connectivity in absence of any structured stimulation. Even after 10000 s of non-modulated external input, the connectivity of excitatory neurons is still slightly modulated according to their PO (Fig. 4D and E right columns).

Figure 4G shows the time evolution of the first Fourier components of the feature connectivity modulation (Fig. 4E) during network growth, stimulation and post-stimulation phase. As previously observed in Fig. 4D and E, there is a drop in average connectivity during the visual stimulation phase, which can be seen in this case as a decrease in the DC component of the modulated signal. Also, during the visual stimulation phase, there is an increase in the first Fourier component, corresponding to the modulation of the signal that happens simultaneously with the decrease in average connectivity (Fig. 4G).

After the stimulation phase, when all excitatory neurons receive once again the same non-modulated external input νext, the average connectivity returns quickly to the value it had before stimulation. The modulation component of the connectivity, on the other hand, decays slowly back to its original value, indicating that there are different time scales for creating and destroying feature specific connectivity based on the modulation of external input. We used non-linear least squares to fit an exponential function to the first Fourier component of the modulated signal during and after stimulation and extracted the time constant for both process, shown in Fig. 4G. As described previously, the process of creating feature specific connectivity is 5 times faster (τ1 ≈ 1000 s) than destroying it (τ2 ≈ 5000 s). The time series of the Fourier components for multiplicity modulation are qualitatively very similar to the feature specific modulation (see Supplementary material).

## Discussion

In our study, we employed a structural plasticity (SP) rule based on the homeostasis of firing rates26,27 to grow random networks of excitatory and inhibitory neurons. We allowed only EE connections to grow, all other synapses were static. We showed that in such a configuration an implementation of homeostatic structural plasticity shares important functional properties with a Hebbian learning rule, with different time scales for the creation and decay of the newly formed associations. In a generic model of visual cortex, we show that feature-specific connectivity, similar to what has been observed in V1 of mice, can emerge from the SP model. This feature-specific connectivity persists even after the feature-specific stimulation has been turned off.

The time scales of structural plasticity in the brain cover several orders of magnitude, ranging from minutes45 up to hours or even days4,9,46. Here we show how the SP model can be adjusted to allow investigations stable through time scales covering at least three orders of magnitude. Minerbi et al.46 continuously imaged cultures of rat cortical neurons and found that, although the distribution of synaptic sizes was stable over days, individual synapses were continuously remodeled. In the SP model, individual pre- and postsynaptic elements are remodeled on a fast time scale, whereas the global network structure evolves at a much slower pace. Therefore, average connectivity and the distributions of in- and outdegrees might be stable, but individual synapses are still plastic, continuously creating new and deleting existing synapses.

We show in Fig. 2 that the time scale of network growth depends on the growth parameter for the synaptic elements (β). A thorough study on the influence of all parameters of the growth rule on network remodeling was beyond the scope of our present study. Therefore, we used the simplest combination of identical linear controllers for both pre- and postsynaptic elements. There is, however, experimental evidence that the growth of axons is somewhat slower than that of spines47, but the exact rules governing the growth of synaptic elements are still unknown. The choice of a linear function to implement the homeostatic controller of firing rates is, in any case, in accordance with empirical studies that demonstrated an increase in the number of newly formed spines on cortical neurons in adult mice after monocular deprivation9 and after small lesions of the retina48. More recently, several studies20,21,22 showed a regulation of firing rate through homeostatic plasticity in the rodent visual cortex in vivo. Keck et al.20 demonstrated an increase of spine size in vivo after a change of sensory input through a retinal lesion, indicating a compensatory recovery through functional plasticity, but no change in spine density after the lesion. This does not completely discard the hypothesis that sensory deprivation could trigger structural plasticity mechanisms, as the homeostatic regulation of activity is probably the result of an interaction between functional and structural plasticity. In this specific case, the spine turnover due to structural plasticity as previously observed by Keck et al.48, and an increase of the strength of existing connections together lead to a recovery of neuronal activity. Regarding presynaptic elements, Canty et al.49 have recently demonstrated axon regrowth in vivo after ablation, with axonal bouton densities similar to the state before the lesion, in accordance to a putative homeostatic mechanism.

Several known aspects of cortical network structure and dynamics were not reproduced in our simulations, such as a broad and skewed distribution of firing rates50,51,52 and synaptic strength53, and a specific motif statistics for pairs and triplets of neurons54,55. The structure and dynamics observed in cortical neuronal networks, however, emerge from the interplay of multiple plasticity mechanisms. A full account of structural remodeling of cortical networks should, therefore, include multiple plasticity mechanisms56,57. With such models, however, is not an easy task to understand what are the effects of individual processes, and what are the effects of combined mechanisms. On one side, by simulating only one plasticity rule in isolation, we were able to report a very interesting property of a growth rule based on firing rate of homeostasis. On the other side, it must remain open how SP interacts with other forms of plasticity, and how the above-mentioned property is changed in the presence of other plastic mechanisms that simultaneously update connectivity in the network.

We have provided support for an idea put forward by Dammasch25: Hebbian plasticity is not necessarily tied to individual synapses, but can also emerge as a system property. The homeostatic control of structural plasticity is achieved on the level of whole neurons, and not of individual synapses. In this model, neurons have control over the number of synaptic elements, which are putative synapses, but not directly over the formation of specific synapses. The realization of a synaptic contact is, in fact, implemented by randomly wiring available elements. Thus, Hebbian learning is implemented through the availability and random wiring of free elements in the network. In contrast to traditional rules implementing Hebbian learning at individual synapses, in the SP model it is not necessary that the neurons keep track of the individual activity of other neurons. Instead, they only need to keep track of their own activity, and a random wiring scheme implements the correlation dependence. Hebbian association, therefore, is formed due to neurons controlling their own total input and output, but not the weight of individual synapses, an idea that may be related to the neurocentric view of learning proposed by Titley et al.58 In their review, Fauth & Tetzlaff59 distinguish two types of structural plasticity rules: (i) Hebbian, if there is an increase (decrease) of the number of synapses during high (low) activity and (ii) homeostatic, if there is an increase (decrease) of the number of synapses during low (high) activity. The SP rule is, according to this definition, a specific variant of homeostatic structural plasticity. On the network level, however, it implements a form of Hebbian plasticity. The classification proposed by Fauth & Tetzlaff59 takes the rules causing the changes in number of pre- and postsynaptic elements into consideration. Another possible classification, however, would consider the effects of the plastic mechanisms on the connectivity. It is actually not an easy task to distinguish these two options in experiments, since what we observe is the effect on the network, and it may be impossible to know what were the mechanisms that led to the observed effects.

Another interesting feature implemented by SP in our specific study are the different time scales for establishing and deleting modulated connectivity. If we consider specific non-random connectivity to implement some sort of memory for previous experiences of the system, this would mean that the system learns faster than it forgets. This appears to happen because the rewiring in the network is triggered by (and depends on) the discrepancy between the actual activity of the neuron and its target rate. Learning takes place if there is modulated external input, and the neuronal firing rate is drawn away from its setpoint, leading to strong rewiring. When the external input is not modulated any more, the neuron’s activity recovers quickly back to its target rate, and the rewiring then becomes very slow, depending on the amplitude of random fluctuations. In a more theoretical framework, Fauth et al.60 showed recently that fast learning and slow forgetting can occur in a stochastic model of structural plasticity. In their model, new synapses are randomly formed with a constant probability, and randomly deleted with a probability that depends on the number of the existing synapses and the current stimulation. In our simulations, in contrast, both the creation and the deletion of synapses depend on the number of synaptic elements of each neuron, which in turn depends on its level of activity. Deletion also depends on the number of existing synaptic contacts between pre- and postsynaptic neurons, due to competition when deleting an existing contact. Different rules for the growth of synaptic elements could of course lead to different dependencies. These rules might also influence other properties of the SP model we describe in this paper, such as the capability to form associations.

Hiratani & Fukai61 demonstrated the formation of cell assemblies of strongly connected cells in a random recurrent network with short-term depression, log-STDP62 and homeostatic plasticity. Similarly to other computational models17,63, the stronger connectivity is accompanied by sustained activity of neurons, consistent with the concept of a working memory. Our results, in contrast, show a memory trace in the connectivity of neurons without sustained activity, more consistent with the idea of contextual memories. Another clear difference between these models concerns the time evolution of average synaptic weights. Hiratani & Fukai61 show an increase in average synaptic weight throughout the stimulation. In our simulations, the average connectivity decreases during stimulation time due to the homeostatic principles underlying the plasticity rule, and it increases after the specific stimulation has been turned off. This would imply that immediately after the end of stimulation, connectivity between the stimulated neurons is lower than baseline. Although this seems counter-intuitive, there have indeed been studies showing perceptual deterioration after trial repetition for subjects tested in a certain task on the same day64,65,66,67, followed by perceptual improvement after 24 and 48 hours67, which would be in agreement to the observed dynamics of connectivity in our simulations. Experimental data on plasticity usually report values of connectivity before and after a stimulation, but do not allow insight into the connectivity dynamics. Knowing the time evolution of these connectivity values during different stimulation protocols could give us important hints about the exact mechanisms of plastic changes and help constrain the plasticity models.

A straight-forward consequence of homeostatic plasticity is the stabilization of activity in neuronal networks. This aspect has been thoroughly studied over many years [see19 for a review on homeostatic plasticity for stabilizing neuronal activity]. The synaptic homeostasis hypothesis should be mentioned in this context68. It states that Hebbian learning during awake states leads to an increase in firing rates, and homeostatic plasticity during sleep states has the goal to restore activity back to baseline levels. Hengen et al.23 recently showed that the opposite is the case. They continuously monitored the firing rate of individual visual cortical neurons in freely behaving rats over several days and showed that homeostasis is actually inhibited by sleep and promoted by wake states. It is of course possible that homeostatic plasticity has an exclusive role for network stabilization, even if it is active during wake and not sleep states. In any case, all these aspects taken together suggest that there could be more to homeostatic plasticity than just stabilizing the network.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## References

1. 1.

Cohan, C. S. & Kater, S. B. Suppression of neurite elongation and growth cone motility by electrical activity. Science 232, 1638–40 (1986).

2. 2.

Van Huizen, F. & Romijn, H. Tetrodotoxin enhances initial neurite outgrowth from fetal rat cerebral cortex cells in vitro. Brain Res. 408, 271–274, https://doi.org/10.1016/0006-8993(87)90386-6 (1987).

3. 3.

Mattson, M. P. & Kater, S. Excitatory and inhibitory neurotransmitters in the generation and degeneration of hippocampal neuroarchitecture. Brain Res. 478, 337–348, https://doi.org/10.1016/0006-8993(89)91514-X (1989).

4. 4.

Harris, K. M. & Kirov, S. A. Dendrites are more spiny on mature hippocampal neurons when synapses are inactivated. Nat. Neurosci. 2, 878–883, https://doi.org/10.1038/13178 (1999).

5. 5.

Wiegert, J. S. & Oertner, T. G. Long-term depression triggers the selective elimination of weakly integrated synapses. Proc. Natl. Acad. Sci. 110, E4510–9, https://doi.org/10.1073/pnas.1315926110 (2013).

6. 6.

Oh, W. C., Parajuli, L. K. & Zito, K. Heterosynaptic Structural Plasticity on Local Dendritic Segments of Hippocampal CA1 Neurons. Cell Rep. 10, 162–169, https://doi.org/10.1016/j.celrep.2014.12.016 (2015).

7. 7.

Zuo, Y., Yang, G., Kwon, E. & Gan, W.-B. Long-term sensory deprivation prevents dendritic spine loss in primary somatosensory cortex. Nature 436, 261–265, https://doi.org/10.1038/nature03715 (2005).

8. 8.

Zuo, Y., Lin, A., Chang, P. & Gan, W.-B. Development of Long-Term Dendritic Spine Stability in Diverse Regions of Cerebral Cortex. Neuron 46, 181–189, https://doi.org/10.1016/j.neuron.2005.04.001 (2005).

9. 9.

Hofer, S. B., Mrsic-Flogel, T. D., Bonhoeffer, T. & Hübener, M. Experience leaves a lasting structural trace in cortical circuits. Nature 457, 313–317, https://doi.org/10.1038/nature07487 (2009).

10. 10.

Xu, T. et al. Rapid formation and selective stabilization of synapses for enduring motor memories. Nature 462, 915–919, https://doi.org/10.1038/nature08389 (2009).

11. 11.

Yang, G., Pan, F. & Gan, W.-B. Stably maintained dendritic spines are associated with lifelong memories. Nature 462, 920–924, https://doi.org/10.1038/nature08577 (2009).

12. 12.

Nikonenko, I., Jourdain, P. & Muller, D. Presynaptic remodeling contributes to activity-dependent synaptogenesis. J. Neurosci. 23, 8498–505 (2003).

13. 13.

Schuemann, A., Klawiter, A., Bonhoeffer, T. & Wierenga, C. J. Structural plasticity of GABAergic axons is regulated by network activity and GABAA receptor activation. Front. Neural Circuits 7, 113, https://doi.org/10.3389/fncir.2013.00113 (2013).

14. 14.

Ninan, I., Liu, S., Rabinowitz, D. & Arancio, O. Early presynaptic changes during plasticity in cultured hippocampal neurons. EMBO J. 25, 4361–4371, https://doi.org/10.1038/sj.emboj.7601318 (2006).

15. 15.

Zheng, P. & Triesch, J. Robust development of synfire chains from multiple plasticity mechanisms. Front. Computat. Neurosci. 8, 66, https://doi.org/10.3389/fncom.2014.00066 (2014).

16. 16.

Miner, D. & Triesch, J. Plasticity-Driven Self-Organization under Topological Constraints Accounts for Non-random Features of Cortical Synaptic Wiring. PLOS Computat. Biol. 12, e1004759, https://doi.org/10.1371/journal.pcbi.1004759 (2016).

17. 17.

Zenke, F., Agnes, E. J. & Gerstner, W. Diverse synaptic plasticity mechanisms orchestrated to form and retrieve memories in spiking neural networks. Nat. Commun. 6, https://doi.org/10.1038/ncomms7922 (2015).

18. 18.

Zenke, F. & Gerstner, W. Hebbian plasticity requires compensatory processes on multiple timescales. Philos. Transact. Royal Soc. B 372, 20160259, https://doi.org/10.1098/rstb.2016.0259 (2017).

19. 19.

Turrigiano, G. Homeostatic synaptic plasticity: local and global mechanisms for stabilizing neuronal function. Cold Spring Harb. perspectives in biology 4, a005736, https://doi.org/10.1101/cshperspect.a005736 (2012).

20. 20.

Keck, T. et al. Synaptic Scaling and Homeostatic Plasticity in the Mouse Visual Cortex In Vivo. Neuron 80, 327–334, https://doi.org/10.1016/j.neuron.2013.08.018 (2013).

21. 21.

Hengen, K. B., Lambo, M. E., Van Hooser, S. D., Katz, D. B. & Turrigiano, G. G. Firing rate homeostasis in visual cortex of freely behaving rodents. Neuron 80, 335–42, https://doi.org/10.1016/j.neuron.2013.08.038 (2013).

22. 22.

Barnes, S. J. et al. Subnetwork-Specific Homeostatic Plasticity in Mouse Visual Cortex In Vivo. Neuron 86, 1290–303, https://doi.org/10.1016/j.neuron.2015.05.010 (2015).

23. 23.

Hengen, K. B., Torrado Pacheco, A., McGregor, J. N., Van Hooser, S. D. & Turrigiano, G. G. Neuronal Firing Rate Homeostasis Is Inhibited by Sleep and Promoted by Wake. Cell 165, 180–191, https://doi.org/10.1016/j.cell.2016.01.046 (2016).

24. 24.

Fox, K. & Stryker, M. Integrating Hebbian and homeostatic plasticity: introduction. Philos. Transactions Royal Soc. B 372, 20160413, https://doi.org/10.1098/rstb.2016.0413 (2017).

25. 25.

Dammasch, I. E. Structural relization of a Hebb-type learning rule. In M, C. R. (ed.) Models of Brain Functions (Cambridge Univeristy Press, 1989).

26. 26.

Butz, M. & van Ooyen, A. A simple rule for dendritic spine and axonal bouton formation can account for cortical reorganization after focal retinal lesions. PLOS Computat. Biol. 9, e1003259, https://doi.org/10.1371/journal.pcbi.1003259 (2013).

27. 27.

Diaz-Pier, S., Naveau, M., Butz-Ostendorf, M. & Morrison, A. Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity. Front. Neuroanat. 10, https://doi.org/10.3389/fnana.2016.00057 (2016).

28. 28.

Ko, H. et al. Functional specificity of local synaptic connections in neocortical networks. Nature 473, 87–91, https://doi.org/10.1038/nature09880 (2011).

29. 29.

Ko, H. et al. The emergence of functional microcircuits in visual cortex. Nature 496, 96–100, https://doi.org/10.1038/nature12015 (2013).

30. 30.

Ko, H., Mrsic-Flogel, T. D. & Hofer, S. B. Emergence of feature-specific connectivity in cortical microcircuits in the absence of visual experience. J. Neurosci. 34, 9812–6, https://doi.org/10.1523/JNEUROSCI.0875-14.2014 (2014).

31. 31.

Demas, J., Eglen, S. J. & Wong, R. O. L. Developmental loss of synchronous spontaneous activity in the mouse retina is independent of visual experience. J. Neurosci. 23, 2851–60 (2003).

32. 32.

Sadeh, S., Clopath, C. & Rotter, S. Emergence of Functional Specificity in Balanced Networks with Synaptic Plasticity. PLOS Computat. Biol. 11, e1004307, https://doi.org/10.1371/journal.pcbi.1004307 (2015).

33. 33.

Gewaltig, M.-O. & Diesmann, M. NEST (NEural Simulation Tool). Scholarpedia 2, 1430, https://doi.org/10.4249/scholarpedia.1430 (2007).

34. 34.

Bos, H. et al. NEST 2.10.0, https://doi.org/10.5281/ZENODO.44222 (2015).

35. 35.

Dammasch, I. E., Wagner, G. P. & Wolff, J. R. Self-stabilization of neuronal networks. Biol. Cybern. 54, 211–222 (1986).

36. 36.

van Ooyen, A. & van Pelt, J. Activity-dependent Outgrowth of Neurons and Overshoot Phenomena in Developing Neural Networks. J. Theor. Biol. 167, 27–43, https://doi.org/10.1006/jtbi.1994.1047 (1994).

37. 37.

Van Ooyen, A., Van Pelt, J. & Corner, M. A. Implications of activity dependent neurite outgrowth for neuronal morphology and network development. J. Theoret. Biol. 172, 63–82, https://doi.org/10.1006/jtbi.1995.0005 (1995).

38. 38.

Butz, M., Van Ooyen, A. & Wörgötter, F. A model for cortical rewiring following deafferentation and focal stroke. Front. Computat. Neurosci. 3, 10, https://doi.org/10.3389/neuro.10.010.2009 (2009).

39. 39.

Butz, M., Steenbuck, I. D. & van Ooyen, A. Homeostatic structural plasticity can account for topology changes following deafferentation and focal stroke. Front. Neuroanat. 8, 115, https://doi.org/10.3389/fnana.2014.00115 (2014).

40. 40.

Butz, M., Steenbuck, I. D. & van Ooyen, A. Homeostatic structural plasticity increases the efficiency of small-world networks. Front. Synaptic Neurosci. 6, 7, https://doi.org/10.3389/fnsyn.2014.00007 (2014).

41. 41.

Tetzlaff, C., Okujeni, S., Egert, U., Wörgötter, F. & Butz, M. Self-Organized Criticality in Developing Neuronal Networks. PLOS Comput. Biol. 6, e1001013, https://doi.org/10.1371/journal.pcbi.1001013 (2010).

42. 42.

Butz, M., Lehmann, K., Dammasch, I. E. & Teuchert-Noodt, G. A theoretical network model to analyse neurogenesis and synaptogenesis in the dentate gyrus. Neural Networks 19, 1490–1505, https://doi.org/10.1016/j.neunet.2006.07.007 (2006).

43. 43.

Butz, M., Teuchert-Noodt, G., Grafen, K. & van Ooyen, A. Inverse relationship between adult hippocampal cell proliferation and synaptic rewiring in the dentate gyrus. Hippocampus 18, 879–898, https://doi.org/10.1002/hipo.20445 (2008).

44. 44.

Brunel, N. Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons. J Comput. Neurosci. 8, 183–208, https://doi.org/10.1023/A:1008925309027 (2000).

45. 45.

Bonhoeffer, T. & Engert, F. Dendritic spine changes associated with hippocampal long-term synaptic plasticity. Nature 399, 66–70, https://doi.org/10.1038/19978 (1999).

46. 46.

Minerbi, A. et al. Long-Term Relationships between Synaptic Tenacity, Synaptic Remodeling, and Network Activity. PLOS Biol. 7, e1000136, https://doi.org/10.1371/journal.pbio.1000136 (2009).

47. 47.

Majewska, A. K., Newton, J. R. & Sur, M. Remodeling of Synaptic Structure in Sensory Cortical Areas In Vivo. J. Neurosci. 26, 3021–3029, https://doi.org/10.1523/JNEUROSCI.4454-05.2006 (2006).

48. 48.

Keck, T. et al. Massive restructuring of neuronal circuits during functional reorganization of adult visual cortex. Nat. Neurosci. 11, 1162–1167, https://doi.org/10.1038/nn.2181 (2008).

49. 49.

Canty, A. J. et al. In-vivo single neuron axotomy triggers axon regeneration to restore synaptic density in specific cortical circuits. Nat. Commun. 4, 777–791, https://doi.org/10.1038/ncomms3038 (2013).

50. 50.

Hromádka, T., DeWeese, M. R. & Zador, A. M. Sparse Representation of Sounds in the Unanesthetized Auditory Cortex. PLOS Biol. 6, e16, https://doi.org/10.1371/journal.pbio.0060016 (2008).

51. 51.

O’Connor, D. H., Peron, S. P., Huber, D. & Svoboda, K. Neural Activity in Barrel Cortex Underlying Vibrissa-Based Object Localization in Mice. Neuron 67, 1048–1061, https://doi.org/10.1016/j.neuron.2010.08.026 (2010).

52. 52.

Mizuseki, K. & Buzsáki, G. Preconfigured, Skewed Distribution of Firing Rates in the Hippocampus and Entorhinal Cortex. Cell Rep. 4, 1010–1021, https://doi.org/10.1016/j.celrep.2013.07.039 (2013).

53. 53.

Cossell, L. et al. Functional organization of excitatory synaptic strength in primary visual cortex. Nature 518, 399–403, https://doi.org/10.1038/nature14182 (2015).

54. 54.

Song, S., Sjöström, P. J., Reigl, M., Nelson, S. & Chklovskii, D. B. Highly nonrandom features of synaptic connectivity in local cortical circuits. PLOS Biol. 3, e68, https://doi.org/10.1371/journal.pbio.0030068 (2005).

55. 55.

Perin, R., Berger, T. K. & Markram, H. A synaptic organizing principle for cortical neuronal groups. Proc. Natl. Acad. Sci. 108, 5419–24, https://doi.org/10.1073/pnas.1016051108 (2011).

56. 56.

Toyoizumi, T., Kaneko, M., Stryker, M. & Miller, K. Modeling the Dynamic Interaction of Hebbian and Homeostatic Plasticity. Neuron 84, 497–510, https://doi.org/10.1016/j.neuron.2014.09.036 (2014).

57. 57.

Deger, M., Seeholzer, A. & Gerstner, W. Multicontact Co-operativity in Spike-Timing?Dependent Structural Plasticity Stabilizes Networks. Cereb. Cortex 61, 247–258, https://doi.org/10.1093/cercor/bhx339 (2017).

58. 58.

Titley, H. K., Brunel, N. & Hansel, C. Toward a Neurocentric View of Learning. Neuron 95, 19–32, https://doi.org/10.1016/j.neuron.2017.05.021 (2017).

59. 59.

Fauth, M. & Tetzlaff, C. Opposing Effects of Neuronal Activity on Structural Plasticity. Front. Neuroanat. 10, 75, https://doi.org/10.3389/fnana.2016.00075 (2016).

60. 60.

Fauth, M., Wörgötter, F. & Tetzlaff, C. Formation and Maintenance of Robust Long-Term Information Storage in the Presence of Synaptic Turnover. PLOS Comput. Biol. 11, e1004684, https://doi.org/10.1371/journal.pcbi.1004684 (2015).

61. 61.

Hiratani, N. & Fukai, T. Interplay between short- and long-term plasticity in cell-assembly formation. PloS one 9, e101535, https://doi.org/10.1371/journal.pone.0101535 (2014).

62. 62.

Gilson, M. & Fukai, T. Stability versus neuronal specialization for STDP: long-tail weight distributions solve the dilemma. PloS one 6, e25339, https://doi.org/10.1371/journal.pone.0025339 (2011).

63. 63.

Sadeh, S., Clopath, C. & Rotter, S. Processing of Feature Selectivity in Cortical Networks with Specific Connectivity. PloS one 10, e0127547, https://doi.org/10.1371/journal.pone.0127547 (2015).

64. 64.

Mednick, S. C., Arman, A. C. & Boynton, G. M. The time course and specificity of perceptual deterioration. Proc. Natl. Acad. Sci. 102, 3881–3885, https://doi.org/10.1073/pnas.0407866102 (2005).

65. 65.

Censor, N., Karni, A. & Sagi, D. A link between perceptual learning, adaptation and sleep. Vision Res. 46, 4071–4074, https://doi.org/10.1016/j.visres.2006.07.022 (2006).

66. 66.

Ofen, N., Moran, A. & Sagi, D. Effects of trial repetition in texture discrimination. Vision Res. 47, 1094–1102, https://doi.org/10.1016/j.visres.2007.01.023 (2007).

67. 67.

Mednick, S., Nakayama, K. & Stickgold, R. Sleep-dependent learning: a nap is as good as a night. Nat. Neurosci. 6, 697–698, https://doi.org/10.1038/nn1078 (2003).

68. 68.

Tononi, G. & Cirelli, C. Sleep and the Price of Plasticity: From Synaptic and Cellular Homeostasis to Memory Consolidation and Integration. Neuron 81, 12–34, https://doi.org/10.1016/j.neuron.2013.12.025 (2014).

## Acknowledgements

Supported by Erasmus Mundus/EuroSPIN, BMBF (grant BFNT 01GQ0830) and DFG (grant EXC 1086). The HPC facilities are funded by the state of Baden-Württemberg through bwHPC and DFG grant INST 39/963-1 FUGG. We thank Sandra Diaz-Pier and Mikaël Naveau from the Research Center Jülich for support on new features of NEST, and Uwe Grauer from the Bernstein Center Freiburg as well as Bernd Wiebelt and Michael Janczyk from the Freiburg University Computing Center for their assistance with HPC applications. The article processing charge was covered by the open access publication fund of the University of Freiburg.

## Author information

### Affiliations

1. #### Bernstein Center Freiburg & Faculty of Biology, University of Freiburg, Freiburg im Breisgau, Germany

• Júlia V. Gallinaro
•  & Stefan Rotter

### Contributions

J.G. and S.R. conceptualized the main goals, J.G. performed the simulations, J.G. and S.R. analyzed the results, S.R. supervised the work, J.G. and S.R. wrote and revised the manuscript.

### Competing Interests

The authors declare no competing interests.

### Corresponding author

Correspondence to Júlia V. Gallinaro.


text:
# Doing Maximum Likelihood Estimation by Hand in R

April 21, 2010
By

(This article was first published on John Myles White: Die Sudelbücher » Statistics, and kindly contributed to R-bloggers)

Lately I’ve been writing maximum likelihood estimation code by hand for some economic models that I’m working with. It’s actually a fairly simple task, so I thought that I would write up the basic approach in case there are readers who haven’t built a generic estimation system before.

First, let’s start with a toy example for which there is a closed-form analytic solution. We’ll ignore that solution and use optimization functions to do the estimation. Starting with this toy example makes it easy to see how well an approximation system can be expected to perform under the best circumstances — and also where it goes wrong if you make poor programming decisions.

Suppose that you’ve got a sequence of values from an unknown Bernoulli variable like so:

 1 2 3  p.parameter <- 0.8 sequence <- rbinom(10, 1, p.parameter) # [1] 0 1 1 1 1 1 1 0 1 0

Given the sequence, we want to estimate the value of the parameter, p, which is not known to us. The maximum likelihood approach says that we should select the parameter that makes the data most probable. For a Bernoulli variable, this is simply a search through the space of values for p (i.e [0, 1]) that makes the data most probable to have observed.

It’s worth pointing out that the analytic solution to the maximum likelihood estimation problem is to use the sample mean. We’ll therefore use mean(sequence) as a measure of the accuracy of our approximation algorithm.

How do we find the parameter numerically? First, we want to define a function that specifies the probability of our entire data set. We assume that each observation in the data is independently and identically distributed, so that the probability of the sequence is the product of the probabilities of each value. For the Bernoulli variables, this becomes the following function:

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  likelihood <- function(sequence, p.parameter) { likelihood <- 1   for (i in 1:length(sequence)) { if (sequence[i] == 1) { likelihood <- likelihood * p.parameter } else { likelihood <- likelihood * (1 - p.parameter) } }   return(likelihood) }

To do maximum likelihood estimation, we therefore only need to use an optimization function to maximize this function. A quick examination of the likelihood function as a function of p makes it clear that any decent optimization algorithm should be able to find the maximum:

 1 2 3 4 5 6 7 8 9 10  possible.p <- seq(0, 1, by = 0.001) jpeg('Likelihood_Concavity.jpg') library('ggplot2') qplot(possible.p, sapply(possible.p, function (p) {likelihood(sequence, p)}), geom = 'line', main = 'Likelihood as a Function of P', xlab = 'P', ylab = 'Likelihood') dev.off()

For single variable cases, I find that it’s easiest to use R’s base function optimize to solve the optimization problem:

 1 2 3 4 5 6 7 8 9 10  mle.results <- optimize(function(p) {likelihood(sequence, p)}, interval = c(0, 1), maximum = TRUE)   mle.results # $maximum # [1] 0.6999843 # #$objective # [1] 0.002223566

Here I’ve used an anonymous function that returns the likelihood of our current data given a value of p; I’ve also specified that the values of p must lie in the interval [0, 1] and asked optimize to maximize the result, rather than minimize, which is the default behavior. Examining the output of optimize, we can see that the likelihood of the data set was maximized very near 0.7, the sample mean. This suggests that the optimization approximation can work. It’s worth noting that the objective value is the likelihood of the data set for the specified value of p. The smallness of the objective for large problems can become a major problem. To understand why, it’s worth seeing what happens as the size of the sample grows from 10 to 2500 samples:

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  error.behavior <- data.frame()   for (n in 10:2500) { sequence <- rbinom(n, 1, p.parameter)   likelihood.results <- optimize(function(p) {likelihood(sequence, p)}, interval = c(0, 1), maximum = TRUE)   true.mle <- mean(sequence)   likelihood.error <- true.mle - likelihood.results$maximum error.behavior <- rbind(error.behavior, data.frame(N = n, Error = likelihood.error, Algorithm = 'Likelihood')) } As you can see, our approximation approach works great until our data set grows, and then it falls apart. This is exactly the opposite of what asymptotical statistical theory tells us should be happening, so it’s clear that something is going very wrong. A quick examination of the results from the last pass through our loop makes clear what’s wrong:  1 2 3 4 5 6 7 8 9 10 11 12  sequence <- rbinom(2500, 1, p.parameter) likelihood.results <- optimize(function(p) {likelihood(sequence, p)}, interval = c(0, 1), maximum = TRUE) likelihood.results #$maximum # [1] 0.9999339 # # $objective # [1] 0 The likelihood of our data is numerically indistinguishable from 0 given the precision of my machine’s floating point values. Multiplying thousands of probabilities together is simply not a viable approach without infinite precision. Thankfully, there’s a very simple solution: replace all of the probabilities with their logarithms. Instead of maximizing the likelihood, we maximize the log likelihood, which involves summing rather than multiplying, and therefore stays numerically stable:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  log.likelihood <- function(sequence, p) { log.likelihood <- 0 for (i in 1:length(sequence)) { if (sequence[i] == 1) { log.likelihood <- log.likelihood + log(p) } else { log.likelihood <- log.likelihood + log(1 - p) } } return(log.likelihood) } You can check that this problem is as easily solved numerically as the original problem by graphing the log likelihood:  1 2 3 4 5 6 7 8 9 10  sequence <- c(0, 1, 1, 1, 1, 1, 1, 0, 1, 0) possible.p <- seq(0, 1, by = 0.001) jpeg('Log_Likelihood_Concavity.jpg') qplot(possible.p, sapply(possible.p, function (p) {log.likelihood(sequence, p)}), geom = 'line', main = 'Log Likelihood as a Function of P', xlab = 'P', ylab = 'Log Likelihood') dev.off() And then you can rerun our error diagnostics using both approaches to confirm that the log likelihood approach does not suffer from the same numerical problems:  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31  error.behavior <- data.frame() for (n in 10:2500) { sequence <- rbinom(n, 1, p.parameter) likelihood.results <- optimize(function(p) {likelihood(sequence, p)}, interval = c(0, 1), maximum = TRUE) log.likelihood.results <- optimize(function(p) {log.likelihood(sequence, p)}, interval = c(0, 1), maximum = TRUE) true.mle <- mean(sequence) likelihood.error <- true.mle - likelihood.results$maximum log.likelihood.error <- true.mle - log.likelihood.results\$maximum   error.behavior <- rbind(error.behavior, data.frame(N = n, Error = likelihood.error, Algorithm = 'Likelihood'), data.frame(N = n, Error = log.likelihood.error, Algorithm = 'Log Likelihood')) }   jpeg('Long-Term_Error_Behavior.jpg') ggplot(error.behavior, aes(x = N, y = Error)) + geom_line(aes(group = Algorithm, color = Algorithm)) + opts(title = 'Long-Term Error Behavior of Two Numerical Approaches', xlab = 'Sample Size', ylab = 'Deviation from True MLE') dev.off()

More generally, given any data set and any model, you can — at least in principle — solve the maximum likelihood estimation problem using numerical optimization algorithms. The general algorithm requires that you specify a more general log likelihood function analogous to the R-like pseudocode below:

 1 2 3 4 5 6 7 8 9 10 11  log.likelihood <- function(sequence.as.data.frame, likelihood.function, parameters) { log.likelihood <- 0   for (i in 1:nrow(sequence.as.data.frame)) { log.likelihood <- log.likelihood + log(likelihood.function(sequence.as.data.frame[i,], parameters)) }   return(log.likelihood) }

Then you need to apply multivariable, constrained optimization tools to find your maximum likelihood estimates. This actually turns out to be a hard problem in general, so I’m going to bail out on the topic here.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

format: plain


text:
Definitions

# life

[lahyf]
life, although there is no universal agreement as to a definition of life, its biological manifestations are generally considered to be organization, metabolism, growth, irritability, adaptation, and reproduction. Protozoa perform, in a single cell, the same life functions as those carried on by the complex tissues and organs of humans and other highly developed organisms. The attributes of life are inherent in such minute structures as viruses, bacteria, and genes, just as they are in the whale and the giant sequoia. In seeking an understanding of life, scientists have broken down many barriers that once separated the physical sciences from the biological sciences; a result of the growth of biochemistry, biophysics, and other interrelated fields of study has been a better understanding of the composition and functioning of living tissues of all kinds.

## Characteristics of Life

Organization is found in the basic living unit, the cell, and in the organized groupings of cells into organs and organisms. Metabolism includes the conversion of nonliving material into cellular components (synthesis) and the decomposition of organic matter (catalysis), producing energy. Growth in living matter is an increase in size of all parts, as distinguished from simple addition of material; it results from a higher rate of synthesis than catalysis. Irritability, or response to stimuli, takes many forms, from the contraction of a unicellular organism when touched to complex reactions involving all the senses of higher animals; in plants response is usually much different than in animals but is nonetheless present. Adaptation, the accommodation of a living organism to its present or to a new environment, is fundamental to the process of evolution and is determined by the individual's heredity. The division of one cell to form two new cells is reproduction; usually the term is applied to the production of a new individual (either asexually, from a single parent organism, or sexually, from two differing parent organisms), although strictly speaking it also describes the production of new cells in the process of growth.

## The Basis of Life

Much of the history of biology and of philosophy as related to biology has been marked by a division of thought between vitalistic (or animistic) and mechanistic (or materialistic) concepts. In the most antithetic interpretations of these concepts, the vitalistic school maintains that there is a vital force that distinguishes the living from the nonliving and the mechanistic school holds that there is no essential difference between the animate and inanimate and that all life can be explained by physical and chemical laws. Such diametrically opposed views have actually seldom been held by investigators of either school; elements of both are usually involved. The animistic school, largely predicated on the inexplicability of the basic phenomena of life, has been greatly overshadowed by the accumulating weight of scientific data. As more and more is learned of the minute details of the structure and composition of the substances that make up the cell (to the extent that some have been synthesized chemically), it has become increasingly apparent that living matter is made up of the same (and only those) elements found in inorganic material, except that they are differently organized.

## The Origin of Life

Fundamental religious concepts center around special creation and belief in the infusion of life into inanimate substance by God or another superhuman entity. On the other hand, many scientists have hypothesized that during an early geological period there gradually formed in the atmosphere increasingly complex organic substances composed of available inorganic compounds and water, utilizing ultraviolet rays and electrical discharges as energy sources. At a certain stage they formed a diffuse solution of "nutrient broth." Then in some way they were drawn together and developed the capacity for self-renewal and self-reproduction. In 1953, S. L. Miller synthesized several of the most basic amino acids in a glass flask by introducing an electrical discharge into an atmosphere of water vapor and some simple compounds thought to have been present naturally at the time when life first developed on earth. A more recent theory now widely held is that life originated in a volcanic setting more than 3.5 billion years ago, perhaps in hot deep-sea vents, utilizing a biochemistry based largely on sulfur and iron. The theory that life on earth came in a simple form from another planet has had small currency, although the discovery by Melvin Calvin of molecules resembling genetic material in meteors has given it some force.

## Bibliography

See M. Calvin, Chemical Evolution (1969); E. Borek, The Sculpture of Life (1973); N. D. Newell, Creation and Evolution (1985); S. W. Fox and K. Dose, Molecular Evolution and the Origins of Life (3d ed. 1990); R. Fortey, Life (1998).

Alcoholic beverage obtained by distillation from wine or other fermented fruit juice or from various cereal grains that have first been brewed. The essential ingredient is usually a natural sugar or a starchy substance that may be easily converted into a sugar. The distillation process is based on the different boiling points of water (212 °F [100 °C]) and alcohol (173 °F [78.5 °C]). The alcohol vapours that arise while the fermented liquid boils are trapped and recondensed to create a liquid of much greater alcoholic strength. The resultant distillate is matured, often for several years, before it is packaged and sold. Seealso aquavit; brandy; gin; liqueur; rum; vodka; and whiskey.

(Latin: “tree of life”) Any of six species (genus Thuja) of resinous, evergreen ornamental and timber conifers of the cypress family, native to North America and eastern Asia. Arborvitae trees or shrubs have thin, scaling outer bark and fibrous inner bark; horizontal or ascending branches; and flattened, spraylike branchlets with scalelike leaves. The Oriental, or Chinese, arborvitae (T. orientalis), a popular ornamental native to Asia, is a gracefully symmetrical shrub. Arborvitae wood is soft and lightweight but very durable, fragrant, and easily worked. The giant arborvitae (T. plicata) is the most important timber-producing species, but the wood of the American arborvitae (T. occidentalis) is also frequently used. False arborvitae (Thujopsis dolabrata) is closely related.

Depiction of inanimate objects for the sake of their qualities of form, colour, texture, composition, and sometimes allegorical or symbolical significance. Still lifes were painted in ancient Greece and Rome. In the Middle Ages they occur in the borders of illuminated manuscripts. The modern still life emerged as an independent genre in the Renaissance. Netherlandish still lifes often depicted skulls, candles, and hourglasses as allegories of mortality, or flowers and fruits to symbolize nature's cycle. Several factors contributed to the rise of still life in the 16th–17th century: an interest in realistic representation, the rise of a wealthy middle class that wanted artworks to decorate its homes, and increased demand for paintings of secular subjects other than portraits in the wake of the Reformation. Dutch and Flemish painters were the masters of still life in the 17th century. From the 18th century until the rise of nonobjective painting after World War II, France was the centre of still-life painting.

In radioactivity, the average lifetime of all the nuclei of a particular unstable atomic species. This time interval is the sum of the lifetimes of all the individual unstable nuclei in a sample, divided by the total number of unstable nuclei present. It is the reciprocal of the decay constant. For a given isotope, the mean life is always 1.443 times its half-life. For example, lead-209 decays to bismuth-209 with a half-life of 3.25 hours and a mean life of 4.69 hours.

Time between birth and death. It ranges from a mayfly's day to certain trees' thousands of years. Its limit appears to depend on heredity, but such factors as (in humans) disease, natural disasters, war, diet, and habits such as smoking reduce it. Maximum life span is theoretical; more meaningful is average life span, which life-insurance companies and actuaries analyze and tabulate. Long-lived progenitors tend to beget long-lived descendants. A very-low-calorie diet appears to prolong life. Reduced infant mortality and improved sanitation and nutrition account for much of the increase since circa 1800—from about 35 to over 70 years in most industrialized countries. The oldest well-documented age reached by a human is 122 years.

Method by which large groups of individuals equalize the burden of financial loss from death by distributing funds to the beneficiaries of those who die. Life insurance is most developed in wealthy countries, where it has become a major channel of saving and investing. There are three basic types of life-insurance contract. Term insurance is issued for a specified number of years; protection expires at the end of the period and there is no cash value remaining. Whole-life contracts run for the whole of the insured's life and also accumulate a cash value, which is paid when the contract matures or is surrendered; the cash value is less than the policy's face value. Endowment contracts run for a specified time period and pay their full face value at the end of the period.

State characterized by the ability to metabolize nutrients (process materials for energy and tissue building), grow, reproduce, and respond and adapt to environmental stimuli. Fossil evidence suggests that earth's first living organisms, bacteria and cyanobacteria, arose about 3.5 billion years ago. All known life-forms possess either DNA or RNA. Viruses, which possess DNA and RNA, cannot reproduce without a host cell and do not metabolize nutrients, and it is uncertain whether they should be classified as living or nonliving. Scientists disagree on the likelihood of extraterrestrial life. Seealso Drake equation.

Interval of time required for one-half of the atomic nuclei of a radioactive sample to decay (change spontaneously into other nuclear species by emitting particles and energy), or the time required for the number of disintegrations per second of a radioactive material to decrease by one-half. Half-lives are characteristic properties of the various unstable atomic nuclei and the particular way in which they decay. Alpha decay and beta decay are generally slower processes than gamma decay.

The half-life of a quantity whose value decreases with time is the interval required for the quantity to decay to half of its initial value. The concept originated in describing how long it takes atoms to undergo radioactive decay, but also applies in a wide variety of other situations.

The term "half-life" dates to 1907. The original term was "half-life period", but that was shortened to "half-life" starting in the early 1950s.

Half-lives are very often used to describe quantities undergoing exponential decay—for example radioactive decay. However, a half-life can also be defined for non-exponential decay processes. For a general introduction and description of exponential decay, see the article exponential decay. For a general introduction and description of non-exponential decay, see the article rate law.

Number of
half-lives
elapsed
Fraction
remaining
Percentage
remaining
0 1/1 100
1 1/2 50
2 1/4 25
3 1/8 12 .5
4 1/16 6 .25
5 1/32 3 .125
6 1/64 1 .563
7 1/128 0 .781
... ... ...
$n$ $1/2^n$ $100\left(1/2^n\right)$
The table at right shows the reduction of the quantity in terms of the number of half-lives elapsed.

## Probabilistic nature of half-life

A half-life often describes the decay of discrete entities, such as radioactive atoms. In that case, it does not work to use the definition "half-life is the time required for exactly half of the entities to decay". For example, if there is just one radioactive atom with a half-life of 1 second, there will not be "half of an atom" left after 1 second. There will be either zero atoms left or one atom left, depending on whether or not the atom happens to decay.

Instead, the half-life is defined in terms of probability. It is the time when the expected value of the number of entities that have decayed is equal to half the original number. For example, one can start with a single radioactive atom, wait its half-life, and measure whether or not it decays in that period of time. Perhaps it will and perhaps it will not. But if this experiment is repeated again and again, it will be seen that it decays within the half life 50% of the time.

In some experiments (such as the synthesis of a superheavy element), there is in fact only one radioactive atom produced at a time, with its lifetime individually measured. In this case, statistical analysis is required to infer the half-life. In other cases, a very large number of identical radioactive atoms decay in the time-range measured. In this case, the central limit theorem ensures that the number of atoms that actually decay is essentially equal to the number of atoms that are expected to decay. In other words, with a large enough number of decaying atoms, the probabilistic aspects of the process can be ignored.

There are various simple exercises that demonstrate probabilistic decay, for example involving flipping coins or running a computer program. See the following websites: , ,

## Formulae for half-life in exponential decay

An exponential decay process can be described by any of the following three equivalent formulae:

$N_t = N_0 \left(1/2\right)^\left\{t/t_\left\{1/2\right\}\right\}$
$N_t = N_0 e^\left\{-t/tau\right\}$
$N_t = N_0 e^\left\{-lambda t\right\}$
where
*$N_0$ is the initial quantity of the thing that will decay (this quantity may be measured in grams, moles, number of atoms, etc.),
*$N_t$ is the quantity that still remains and has not yet decayed after a time t,
*$t_\left\{1/2\right\}$ is the half-life of the decaying quantity,
*τ is a positive number called the mean lifetime of the decaying quantity,
*λ is a positive number called the decay constant of the decaying quantity.
The three parameters $t_\left\{1/2\right\}$, τ, and λ are all directly related in the following way:
$t_\left\{1/2\right\} = frac\left\{ln \left(2\right)\right\}\left\{lambda\right\} = tau ln\left(2\right)$
where ln(2) is the natural logarithm of 2 (approximately 0.693).

By plugging in and manipulating these relationships, we get all of the following equivalent descriptions of exponential decay, in terms of the half-life:

$N_t = N_0 \left(1/2\right)^\left\{t/t_\left\{1/2\right\}\right\} = N_0 2^\left\{-t/t_\left\{1/2\right\}\right\} = N_0 e^\left\{-tln\left(2\right)/t_\left\{1/2\right\}\right\}$
$t_\left\{1/2\right\} = t/log_2\left(N_0/N_t\right) = t/\left(log_2\left(N_0\right)-log_2\left(N_t\right)\right) = \left(log_\left\{2^t\right\}\left(N_0/N_t\right)\right)^\left\{-1\right\} = tln\left(2\right)/ln\left(N_0/N_t\right)$
Regardless of how it's written, we can plug into the formula to get

• $N_t=N_0$ at t=0 (as expected—this is the definition of "initial quantity")
• $N_t=\left(1/2\right)N_0$ at $t=t_\left\{1/2\right\}$ (as expected—this is the definition of half-life)
• $N_t$ approaches zero when t approaches infinity (as expected—the longer we wait, the less remains).

### Decay by two or more processes

Some quantities decay by two exponential-decay processes simultaneously. In this case, the actual half-life T1/2 can be related to the half-lives t1 and t2 that the quantity would have if each of the decay processes acted in isolation:
$frac\left\{1\right\}\left\{T_\left\{1/2\right\}\right\} = frac\left\{1\right\}\left\{t_1\right\} + frac\left\{1\right\}\left\{t_2\right\}$
For three or more processes, the analogous formula is:
$frac\left\{1\right\}\left\{T_\left\{1/2\right\}\right\} = frac\left\{1\right\}\left\{t_1\right\} + frac\left\{1\right\}\left\{t_2\right\} + frac\left\{1\right\}\left\{t_3\right\} + cdots$
For a proof of these formulae, see Decay by two or more processes.

### Examples

There is a half-life describing any exponential-decay process. For example:

• The current flowing through an RC circuit or RL circuit decays with a half-life of $RCln\left(2\right)$ or $ln\left(2\right)L/R$, respectively.
• In a first-order chemical reaction, the half-life of the reactant is $ln\left(2\right)/lambda$, where λ is the reaction rate constant.
• In radioactive decay, the half-life is the length of time after which there is a 50% chance that an atom will have undergone nuclear decay. It varies depending on the atom type and isotope, and is usually determined experimentally.

## Half-life in non-exponential decay

Many quantities decay in a way not described by exponential decay—for example, the evaporation of water from a puddle, or (often) the chemical reaction of a molecule. In this case, the half-life is defined the same way as before: The time elapsed before half of the original quantity has decayed. However, unlike in an exponential decay, the half-life depends on the initial quantity, and changes over time as the quantity decays.

As an example, the radioactive decay of carbon-14 is exponential with a half-life of 5730 years. If you have a quantity of carbon-14, half of it (on average) will have decayed after 5730 years, regardless of how big or small the original quantity was. If you wait another 5730 years, one-quarter of the original will remain. On the other hand, the time it will take a puddle to half-evaporate depends on how deep the puddle is. Perhaps a puddle of a certain size will evaporate down to half its original volume in one day. But if you wait a second day, there is no reason to expect that precisely one-quarter of the puddle will remain; in fact, it will probably be much less than that. This is an example where the half-life reduces as time goes on. (In other non-exponential decays, it can increase instead.)

For specific, quantitative examples of half-lives in non-exponential decays, see the article Rate equation.

A biological half-life is also a type of half-life associated with a non-exponential decay, namely the decay of the activity of a drug or other substance after it is introduced into the body.