Homeostatic structural plasticity increases the efficiency of small-world networks

In networks with small-world topology, which are characterized by a high clustering coefficient and a short characteristic path length, information can be transmitted efficiently and at relatively low costs. The brain is composed of small-world networks, and evolution may have optimized brain connectivity for efficient information processing. Despite many studies on the impact of topology on information processing in neuronal networks, little is known about the development of network topology and the emergence of efficient small-world networks. We investigated how a simple growth process that favors short-range connections over long-range connections in combination with a synapse formation rule that generates homeostasis in post-synaptic firing rates shapes neuronal network topology. Interestingly, we found that small-world networks benefited from homeostasis by an increase in efficiency, defined as the averaged inverse of the shortest paths through the network. Efficiency particularly increased as small-world networks approached the desired level of electrical activity. Ultimately, homeostatic small-world networks became almost as efficient as random networks. The increase in efficiency was caused by the emergent property of the homeostatic growth process that neurons started forming more long-range connections, albeit at a low rate, when their electrical activity was close to the homeostatic set-point. Although global network topology continued to change when neuronal activities were around the homeostatic equilibrium, the small-world property of the network was maintained over the entire course of development. Our results may help understand how complex systems such as the brain could set up an efficient network topology in a self-organizing manner. Insights from our work may also lead to novel techniques for constructing large-scale neuronal networks by self-organization.

1. Introduction

The synaptic wiring of cortical networks is key to the functionality of the brain and a precondition for all cognitive behavior (Park and Friston, 2013). How are synaptic connections set up between the brain's billions of neurons so that cost efficiency (Latora and Marchiori, 2001), for example in terms of wiring length or energy consumption associated with information transmission, is maximized? Due to the presence of long-range connectivity, the brain can be regarded as a highly efficient graph, in the sense that information has to pass only a few intermediate neurons (“nodes”) to travel across the whole brain (Kaiser and Hilgetag, 2006). At the same time, brain connectivity is local and clustered (Hilgetag and Kaiser, 2004), making networks less vulnerable because of the many alternative routes that exist between two nodes. Network topology that combines these two properties, long-range connectivity and high clustering, is called small-world (Watts and Strogatz, 1998). Small-world networks have the advantages of local connectivity combined with a high efficiency brought about by a small number of long-range connections (Watts and Strogatz, 1998). The brain seems to be optimized for maximizing cost efficiency of parallel information processing (Achard and Bullmore, 2007) by widely adopting small-world topology (Sporns and Zwi, 2004; Bassett and Bullmore, 2006). Already in the infant brain, connectivity has small-world characteristics (Fransson et al., 2011). A recent fMRI study reported an increase in small-worldness of brain networks in the first 2 years of life that goes along with a growing number of long-distance connections and therefore an increase in global efficiency (Gao et al., 2011). So far, analysis of topology has focused on neuronal networks with static connectivity, in which plasticity arises from changes in the strength of existing synapses rather than from the rewiring of connectivity. However, particularly during development but also in adulthood, connectivity is not fixed (Butz et al., 2009b), and synapse formation goes along with massive synapse deletion and reorganization of connectivity (Missler et al., 1993; Siegel and Lohmann, 2013). This observation triggered the question of this study: how does network topology develop during ontogeny when synapse formation and pruning cause a constant rewiring of network connectivity?

Various explanations have been proposed to account for the development of synaptic connectivity. For example, axon guidance molecules may form the basis of a genetically-encoded developmental scheme (Yamamoto et al., 2002; Borisyuk et al., 2008). Target neurons may secret signaling molecules that can attract or repel axons. Axons can then follow or move away from the concentration gradient. This form of chemotaxis is usually discussed in the context of the formation of global connectivity. However, for the formation of local connectivity chemical cues are less suitable, since they fail to establish stable gradients over very short distances, below 0.7 mm (Kaiser et al., 2009). Synaptic adhesion molecules (e.g., neuroligins) were proposed as molecular cues for local synapse formation (Scheiffele et al., 2000; Stan et al., 2010). In addition, mechanical forces in the tissue could influence neurite outgrowth (Franze, 2013). Alternatively, the formation of local connectivity may basically be random (Braitenberg and Schüz, 1998), just depending on the accidental overlap of axons and dendrites (Binzegger et al., 2004; van Pelt and van Ooyen, 2013; McAssey et al., 2014; van Ooyen et al., 2014). With random synapse formation, the chance of forming connections decreases with increasing distance between neurons (Kaiser et al., 2009).

Although random overlap of axons and dendrites may explain emerging connectivity, it does not account for the actual driving forces underlying neurite outgrowth. There is ample evidence that electrical activity shapes neuronal morphology (Dalva et al., 1994; Wong and Ghosh, 2002; Uesaka et al., 2005; Butz et al., 2009b) and network formation (Ko et al., 2013). Electrical activity influences neurite outgrowth and retraction (McKinney et al., 1999a; Konur and Ghosh, 2005; Lohmann and Wong, 2005; Tailby et al., 2005; Hutchins and Kalil, 2008), as well as the formation and deletion of axonal boutons and dendritic spines (McKinney et al., 1999b; Groc et al., 2002; Jourdain et al., 2003; Kirov et al., 2004; Hofer et al., 2009). Experimental findings further suggest that neuronal morphogenesis is driven by the need of neurons to establish and maintain a homeostatic equilibrium of their average electrical activity (Kirov et al., 2004; Keck et al., 2008, 2013). Restoration of neuronal firing rate after a change in neuronal input has been found experimentally after, for example, focal retinal lesions (Hengen et al., 2013). Based on these experimental findings, we postulated that whenever during development (van Ooyen and van Pelt, 1994; Van Ooyen et al., 1995; Tetzlaff et al., 2010) or in the mature brain (Butz et al., 2008; Butz and van Ooyen, 2013) a neuron senses a deviation of its electrical activity from a homeostatic set-point, it will initiate changes in its morphology that increase the chance of synapse formation or break existing connections so that its firing rate may be restored. Here we investigate what the impact is of neurons regulating their electrical activities homeostatically on network formation and emerging network topologies.

In order to study the impact of electrical activity on emerging network topology, we used our recent Model of Structural Plasticity (MSP) (Butz and van Ooyen, 2013). There are important earlier models of homeostatic structural plasticity, such as the compensation model by Dammasch et al. (Dammasch et al., 1986, 1988; Cromme and Dammasch, 1989; Butz and Teuchert-Noodt, 2006; Butz et al., 2006) and the activity-dependent neurite outgrowth model by van Ooyen (van Ooyen and van Pelt, 1994; Van Ooyen et al., 1995). The latter model, which studied the reciprocal interactions between neuronal activity and network formation, successfully accounted for experimental data on developing cell cultures (van Ooyen and van Pelt, 1994; van Oss and van Ooyen, 1997; Abbott and Rohrkemper, 2007; Tetzlaff et al., 2010). However, both models are not suited for studying topology development. The compensation model lacks topology at all, whereas the representation of neuritic fields by circles in the model by van Ooyen imposes too strong constraints on network topology. In van Ooyen's model, neurons always connect to their direct neighbors before connecting to more distant neurons. Therefore, we developed our MSP (Butz and van Ooyen, 2013), in which we replaced the circle representation by discrete synaptic elements whose numbers change in an activity-dependent manner. Synapses are formed in a random and distance-dependent way by combining synaptic elements from different neurons.

In MSP, the local activity-dependent growth process in combination with a simple kernel function favoring the formation of short-range connections over long-range connections shapes the development of small-world networks. Our simulation results revealed an interesting property of homeostatic growth: as soon as most neurons approached homeostasis in electrical activity, they started forming more long-range connections than expected from the kernel function. Although the clustering coefficient decreased as a result of the formation of long-range connections, the network maintained its small-world property. Furthermore, connectivity became more diverse, as indicated by a decreasing betweenness centrality, and attained a higher global efficiency (defined as the averaged inverse of the shortest paths between all neurons in the network) than small-world networks without homeostasis. Our findings may account for experimental data on the topology of developing dissociated cell cultures (Downes et al., 2012). Interesting similarities were also found between the model and early human brain development (Gao et al., 2011) with respect to an increasing number of long-range connections and an increasing global efficiency during development.

2. Materials and Methods

2.1. Neuron Model

The model network consisted of n = 400 neurons, of which 80% were excitatory and 20% inhibitory. Excitatory neurons were placed with some jitter on a 20 × 16 grid with a spatial distance between two grid points of 150 μm. Inhibitory neurons were placed evenly between the excitatory neurons. The same network layout was used as in Butz and van Ooyen (2013). We used Izhikevich's model (Izhikevich, 2003) to simulate neuronal electrical activity. This model has two differential equations, one for the membrane potential v (in mV) and one for a recovery variable u (in mVms⁻¹), enabling re-polarization after an action potential:

$\begin{array}{cc}\begin{array}{l}\frac{dv}{dt}=k_1{v}^2+k_{2}v+k_3-u+I\\ \frac{du}{dt}=a(bv-u)\hfill \end{array}& \left(1\right)\end{array}$

where k₁ = 0.04 mV⁻¹ ms⁻¹, k₂ = 5 ms⁻¹, k₃ = 140 mVms⁻¹, and t is in ms. Every time a neuron fires (v ≥ 30 mV), v and u are reset:

$\begin{array}{cc}\text{if}v\ge 30 \text{mV},\text{then}\{\begin{array}{ccc}v\hfill & \leftarrow \hfill & c\hfill \\ u\hfill & \leftarrow \hfill & u+d\hfill \end{array}& \left(2\right)\end{array}$

As in Izhikevich (2003), the following parameter values were used: a = 0.1 ms⁻¹, b = 0.2 ms⁻¹, c = −65 mV, d = 2 mVms⁻¹. We used the same dynamics for excitatory and inhibitory neurons. Each neuron receives input I = I_syn + I_ext, consisting of synaptic input I_syn from within the network and external input I_ext. Neurons interchange electrical signals on a millisecond timescale without synaptic delay. Synaptic input consists of the incoming action potentials from the presynaptic neurons low-pass filtered by an exponential filter function $h(t)=exp\left(-\frac{t}{\mu }\right)$ with decay constant μ = 5 ms. Network connectivity W_i,j is defined as the number of synapses from neuron j to i. If a synapse exists, it has a fixed strength of 1 mVms⁻¹. Neurons are either excitatory or inhibitory. Indices refer to excitatory neurons if i or j ∈ {Ex} and to inhibitory neurons if i or j ∈ {In}. As in Izhikevich (2003) and Butz and van Ooyen (2013), I_ext is delivered as white noise with mean 5 mVms⁻¹ and standard deviation 1 mVms⁻¹. In some cases lower input values were required and specified where relevant.

The intracellular calcium concentration of a neuron is used as a low-passed filtered average of its firing frequency (Butz and van Ooyen, 2013). Every time a neuron fires, the calcium concentration is increased by a fixed amount; otherwise the concentration decreases exponentially to zero:

$\begin{array}{cc}\frac{d{\left[C{a}^{2+}\right]}_i}{dt}=\{\begin{array}{cc}-\frac{{\left[C{a}^{2+}\right]}_i}{{\tau }_{Ca}}+\beta \hfill & \text{if}v\ge 30\text{mV}\hfill \\ -\frac{{\left[C{a}^{2+}\right]}_i}{{\tau }_{Ca}}\hfill & \text{otherwise}\hfill \end{array}& \left(3\right)\end{array}$

where β = 0.001 ms⁻¹ and τ_Ca = 10,000 ms. We defined a set-point ϵ = 0.7 in calcium concentration, corresponding to an intermediate level of average electrical activity. Every time the neuron's calcium concentration deviates from ϵ, it will induce structural changes in connectivity to restore the desired average level of activity, as described below.

2.2. Synapse Model for Structural Plasticity

During development, neurons show a pronounced formation and pruning of synapses. To simulate reorganization of synaptic connectivity, standard models of synaptic plasticity, in which connectivity is considered fixed with plasticity merely arising from changes in the strength of existing synapses (modeled as weight factors), are not suitable. For this study on the development of neuronal networks, we therefore used our Model of Structural Plasticity (MSP) (Butz and van Ooyen, 2013). The characteristic feature of this model is that it represents synapses as consisting of two recombinable synaptic elements, namely an axonal element and a dendritic element. Axonal elements represent countable presynaptic specializations for transmitter release such as axonal terminals or boutons, while dendritic elements represent the postsynaptic counterparts, i.e., postsynaptic receptor plates on dendritic spines or the dendritic shaft. Axonal and dendritic elements are either excitatory or inhibitory. Any neuron i can form A_i axonal elements, which are excitatory if the neuron is excitatory, or inhibitory if the neuron is inhibitory. At the same time, a neuron, irrespective of its type, can express D^ex_i excitatory and Dⁱⁿ_i inhibitory dendritic elements. Complementary elements can merge to form a synapse (excitatory axonal with excitatory dendritic elements, and inhibitory axonal with inhibitory dendritic elements). Synaptic elements form and delete independently from a synaptic contact partner. In case a neuron deletes a synaptic element that is bound in a synapse, the complementary synaptic element on the other neuron remains and becomes vacant and available again for synapse formation with a new target. Therefore MSP allows for synaptic rewiring.

2.3. Homeostatic Growth Rules

It is well documented that neurons change their morphology in an activity-dependent fashion during development (Butz et al., 2009b; van Ooyen, 2011). The way in which neurons change their morphology suggests that they try to maintain homeostasis of electrical activity (van Ooyen and van Pelt, 1994; Van Ooyen et al., 1995; Butz et al., 2009a; van Ooyen, 2011; Butz and van Ooyen, 2013). That is, neurons in which the average activity is too low start forming new neuritic structures, whereas neurite growth is halted or neurites are pruned when activity is higher than a desired level (homeostatic set-point). In our MSP, we abstracted away from describing detailed neuronal morphology and assumed that changes in morphology effectively result in changing numbers of axonal and dendritic elements—the contact sites for synaptic connections. Homeostatic adaptation of the number of axonal and dendritic elements was modeled by the following growth rule (Figure 1), in which dz represents the change in the number of synaptic elements, with z being A, D^ex or Dⁱⁿ:

$\begin{array}{cc}\frac{dz}{dt}=\nu \left(\frac{2}{1+e^{\left(\left[C{a}^{2+}\right]-ϵ\right)/0.1}}-1\right)& \left(4\right)\end{array}$

where ϵ is the set-point of the intracellular calcium concentration, corresponding to a desired average firing rate of the neuron, and ν is the growth rate of synaptic elements. We chose ν = 10⁻⁴ ms⁻¹, which is slow enough that electrical dynamics and structural dynamics do not interfere, yet fast enough not to slow down the simulations unnecessarily. For reasons of simplicity and because of a lack of detailed experimental data, we applied identical sigmoid growth rules to all types of synaptic elements. In a random and distance-dependent recombination process, newly formed synaptic elements were distributed to matching synaptic elements, thereby forming synapses and creating the pattern of connectivity in the network.

FIGURE 1

Figure 1. Identical sigmoidal growth rules were used for all types of synaptic elements to determine the change dz/dt in the number of synaptic elements in dependence on the intracellular calcium concentration [Ca²⁺]. z needs to be replaced by the respective type of element A, D^ex or Dⁱⁿ. The homeostatic set-point ϵ is the value of the calcium concentration where dz/dt = 0.

2.4. Kernel Function for Synapse Formation

As in our previous work on MSP (Butz and van Ooyen, 2013), we assumed that synapse formation is more likely between adjacent neurons than between distant ones. We applied a two-dimensional Gaussian-shaped kernel function centered at the x, y-coordinates of neuron i, with K_i,j the distance-dependent likelihood for synapse formation between neuron j and i:

$\begin{array}{cc}{K}_{i,j,i\ne j}=e^{-\frac{{\left(po{s}_{xj}-po{s}_{xi}\right)}^2+{\left(po{s}_{yj}-po{s}_{yi}\right)}^2}{{\sigma }^2}}& \left(5\right)\end{array}$

where pos_xi is the x-coordinate and pos_yi is the y-coordinate of postsynaptic neuron i, and pos_xj and pos_yj are the coordinates of presynaptic neuron j. The probability for autapse connections (i.e., a neuron connecting to itself) was set to zero (K_i,j = 0 for i = j). For these simulations we chose σ = 1 × 150 μm where 150 μm is the distance between two grid points.

In order to investigate the impact of this distance-dependency on emerging network topologies, we additionally grew networks with a flat kernel, i.e., with K_{i,j, i ≠ j} = 1 and K_i,j = 0 for i = j.

2.5. Synapse Formation and Deletion

The MSP algorithm (Butz and van Ooyen, 2013) proceeded in three steps to update network connectivity in an activity-dependent fashion. First, electrical activity of all neurons was computed continuously over time. Secondly, depending on the average electrical activity and according to the homeostatic growth rule (Equation 4), the number of elements changed continuously, too. Thirdly, at discrete time points, network connectivity W was updated by synapse formation and deletion. Because of the low growth rate ν = 10⁻⁴ ms⁻¹, changes in connectivity are very slow compared to changes in activity, so it was not necessary to update connectivity at every time step but only at every 100 ms. The timescale of network formation in the model corresponds in principle to a timescale of days or weeks. However, as described above under homeostatic growth rules, the timescale of structural changes was chosen so that the simulations were not unnecessarily slowed down, which means that the total duration of the simulation in milliseconds is not required to sum up to weeks.

2.5.1. Synapse deletion

Since network connectivity is updated at discrete time steps but synaptic elements change continuously over time due to the activity-dependent growth rules, it can happen that a neuron has more outgoing synapses than axonal elements or more incoming synapses than dendritic elements at the time of the next update in network connectivity. In this case, the neuron has to delete the surplus of synapses and to update connectivity.

To update connectivity, the algorithm needs to select which synapses are to be removed. All synapses have an equal chance of being deleted. Note, however, that multiple synapses can co-exist from neuron j to i and that the more synapses there are, the higher the chance that a synapse between neuron j and i will be deleted. The probability P^del_i,j for synapse deletion between neuron j and i is computed by the following master equation that captures four different cases:

$\begin{array}{cc}{P}_{i,j}^{del}=\frac{W_{i,j}}{{\displaystyle \sum W_{k,l}}}& \left(6\right)\end{array}$

For deletion of incoming synapses, we need to distinguish between excitatory and inhibitory synapses in Equation (6). For deleting incoming excitatory synapses of neuron i ∈ {In ∪ Ex}, we sum up W_k,l over all l ∈ {Ex}. For deleting incoming inhibitory synapses of neuron i ∈ {In ∪ Ex}, we sum up W_k,l over all l ∈ {In}. For deletion of outgoing excitatory synapses of excitatory presynaptic neuron j with j ∈ {Ex}, in Equation (6) all synapses are considered to any postsynaptic neuron k with k ∈ {In ∪ Ex}. Thus, we sum up W_k,l over all k ∈ {In ∪ Ex}. The same holds true for outgoing inhibitory synapses with j ∈ {In}.

Sequentially, outgoing and incoming excitatory and inhibitory synapses were selected for deletion. For every type of synapse, the accumulated sum of P^del_i,j [see description of Equation (6) for the range of i and j] gave a probability distribution from which we drew the required number of synapses to be deleted. The selected synapse was deleted by reducing the respective entry W_i,j in the connectivity matrix by one. It can happen that more then one synapse is selected for deletion from the same connection j to i. In this case, the implementation of the algorithm makes sure that the number of synapses to be deleted did not exceed W_i,j. Whenever a neuron deletes a synaptic element that is bound in a synapse, the complementary synaptic element on the other neuron remains and becomes vacant again.

2.5.2. Synapse formation

For synapse formation, the algorithm checked whether a neuron gained vacant synaptic elements, i.e., whether the total number of synaptic elements exceeded the number of bound synaptic elements of this type. Matching vacant synaptic elements (vacant excitatory axonal elements A^ex,vac_j with vacant excitatory dendritic elements D^ex,vac_i, and inhibitory axonal with inhibitory dendritic elements D^in,vac) were randomly distributed among each other with probability density function P^form. The probability P^form_i,j for forming new synapses between neuron j and i depended on the number of vacant synaptic elements they offered and on the Euclidean distance between neuron j and i:

$\begin{array}{cc}{P}_{i,j}^{form}=\left\{\begin{array}{cc}j\in \left\{Ex\right\}:\hfill & \frac{A_j^{vac}{D}_i^{ex,vac}}{{\displaystyle {\sum }_{\iota \in \left\{Ex\right\}}{A}_{\iota }^{vac}}{\displaystyle {\sum }_{\kappa \in \{Ex\cup In\}}{D}_{\kappa }^{ex,vac}}}{K}_{ij}\hfill \\ j\in \left\{In\right\}:\hfill & \frac{A_j^{vac}{D}_i^{in,vac}}{{\displaystyle {\sum }_{\iota \in \left\{In\right\}}{A}_{\iota }^{vac}}{\displaystyle {\sum }_{\kappa \in \{Ex\cup In\}}{D}_{\kappa }^{in,vac}}}{K}_{ij}\hfill \end{array}\right\}\text{ with}i\in \{Ex\cup In\}.& \left(7\right)\end{array}$

The minor number of vacant excitatory and inhibitory axonal or dendritic elements determined how many new excitatory and inhibitory synapses, respectively, could at most be formed (so-called potential synapses) in every update of connectivity. Thus, the number of excitatory and inhibitory potential synapses equaled

$\begin{array}{cc}\begin{array}{c}{M}^{PotSyn,ex}=min\left({\displaystyle {\sum }_{\iota \in \left\{Ex\right\}}{A}_{\iota }^{vac}},{\displaystyle {\sum }_{\kappa \in \{Ex\cup In\}}{D}_{\kappa }^{ex,vac}}\right)\\ M^{PotSyn,in}=min\left({\displaystyle {\sum }_{\iota \in \left\{In\right\}},}{A}_{\iota }^{vac},{\displaystyle {\sum }_{\kappa \in \{Ex\cup In\}}{D}_{\kappa }^{in,vac}}\right)\end{array}& \left(8\right)\end{array}$

for every update in connectivity.

From this distribution, the algorithm chose at maximum M^PotSyn,ex excitatory and M^PotSyn,in inhibitory connections at which new synapses were created. The respective entries W_i,j in the connectivity matrix were then increased by one. A connection was chosen by drawing a random number from a uniform distribution and comparing it to the accumulated probabilities P^form_i,j for all excitatory connections and all inhibitory connections of the entire network. That connection was chosen that had the highest accumulated probability that the random number just did not exceed. If, for this try, the random number exceeded all accumulated probabilities, no synapse was formed. Hence, not necessarily all of the potential synapses were formed.

Additionally, synapse formation needed to fulfil the condition that the number W⁺_i,j of newly-formed synapses from neuron j to i did not exceed the number of vacant synaptic elements that neuron j and i offered:

$\begin{array}{cc}{W}_{i,j}^{+}\le \left\{\begin{array}{lll}j\in \left\{Ex\right\}\hfill & :\hfill & \min (A_j^{vac},D_i^{ex,vac})\hfill \\ j\in \left\{In\right\}\hfill & :\hfill & \min (A_j^{vac},D_i^{in,vac})\hfill \end{array}\right\}\text{ with}i\in \{Ex\cup In\}.& \left(9\right)\end{array}$

In every update, this condition was checked and synapse formation infringing this condition was rejected. Alternatively, the update of connectivity can also be implemented in a purely local fashion (Butz and van Ooyen, 2013). For small networks, the current implementation is more efficient than the original MSP algorithm (Butz and van Ooyen, 2013) because run time is not dependent on the growing numbers of vacant synaptic elements but proportional to the square of the number of neurons in the network. However, for large networks with n >> 1000, n² will quickly out-range the number of vacant synaptic elements, in which case looping over the number of synaptic elements would be faster. Particularly in Matlab, the current description of MSP allows for an elegant vectorized implementation providing additional speed up.

2.6. Development of Non-Homeostatic Networks

To investigate the impact of homeostasis in electrical activity on developing network topology, we created for every homeostatic network a corresponding non-homeostatic network. At every update in connectivity, we took the number of synapses from the homeostatic network and distributed them in the non-homeostatic network with the same kernel function as used in the homeostatic network. Hence, the placement of synapses in the non-homeostatic network was purely dependent on the kernel function but did not meet the homeostasis criterion. The algorithmic implementation for placing synapses was the same for homeostatic and non-homeostatic networks, with P^form_i,j = K_i,j. Instead of distributing M^PotSyn synapses (Equation 8), we simply distributed the total number of synapses from the homeostatic network. Since synapse formation was not limited by numbers of vacant elements, Equation (9) was not applicable.

2.7. Topology Measurements

A neuronal network can be seen as a graph with the neurons being the nodes and the synapses being the edges or links between two nodes. Since the presynaptic neuron always activates the postsynaptic neuron (and never the other way around), we regard the graph as directed. At every update in connectivity, we assessed those graph theoretical measures that are indicative of small-worldness and network efficiency. In addition, betweenness centrality was measured to determine the importance of nodes in the network. To reduce the complexity of the assessment, we considered only the topology of excitatory synaptic connections W^ex,ex between the n^ex excitatory neurons. For the graph theoretical assessments, the brain connectivity toolbox by Olaf Sporns et al. (Rubinov and Sporns, 2010) was used.

2.7.1. Weighted characteristic path length

The characteristic path length L is the average shortest path from one node to any other node in the network. Path length is defined as the number of links that need to be traveled to go from one node (possibly via intermediate nodes) to any other node. On top of this definition, a direct link between two nodes in a weighted network is considered “shorter” the stronger the weight of the link is. For our network, we take the number of synapses W^ex,ex_i,j between two directly linked neurons j and i, with i, j ∈ {Ex}, as the weight of the connection and the inverse 1/W^ex,ex_i,j as the length l_i,j of the connection. The shortest path d_i,j is then the smallest sum of connection lengths that lead from neuron j to i via any intermediate neurons. Our calculation of the weighted characteristic path length was based on an implementation of Dijkstra's algorithm for computing the shortest path in weighted directed networks by Rubinov and Sporns (2010).

2.7.2. Weighted clustering coefficient

The clustering coefficient is an indication of how strongly nodes in a network are interconnected. It can be measured by the number of triangles, ${\tilde{t}}_i^D$, one node forms with any other two nodes in the network divided by the total number of possible triangles, T^D_i. For weighted and directed graphs, one needs to realize that the adjacency matrix (in our case, the connectivity matrix W^ex,ex of the excitatory neurons) is not symmetric and that the entries of the adjacency matrix are not one but can have any weight. According to Fagiolo (2007), the clustering coefficient of a single node in a weighted directed network ${\tilde{C}}_i^D$ is computed as

$\begin{array}{cc}{\tilde{C}}_i^D\left(W^{ex,ex}\right)=\frac{{\tilde{t}}_i^D}{T_i^D}=\frac{{\left[{\left(W^{ex,ex}\right)}^{\left[1/3\right]}+{\left({\left(W^{ex,ex}\right)}^T\right)}^{\left[1/3\right]}\right]}_{i,i}^3}{2\left[d_i^{tot}\left(d_i^{tot}-1\right)-2d_i^{\leftrightarrow }\right]}& \left(10\right)\end{array}$

where (W^ex,ex)^[1/k] = w^1/k_i,j, the k^th root of the entries of the matrix for i,j ∈ {Ex}, and (W^ex,ex)^T is the transposed W^ex,ex matrix. The variable d^tot_i denotes the total degree of node i (the degree counts the number of either incoming or outgoing edges per node, and the total degree is the sum over both the incoming and outgoing edges), and d^↔_i stands for the number of bilateral edges of node i (the number of nodes node i projects to and receives edges from, excluding self-interactions of node i).

The overall clustering coefficient of the network is thus ${\tilde{C}}^D={(n^{ex})}^{-1}{\displaystyle {\sum }_{i=1}^N{\tilde{C}}_i^D}$. Note that for this assessment we only considered n^ex excitatory nodes and excitatory connections W^ex,ex. For a more detailed description of Equation (10), see Fagiolo (2007). We computed the clustering coefficient of the developing neuronal networks at every update in connectivity by using the implementation by Mika Rubinov from the brain connectivity toolbox (Rubinov and Sporns, 2010).

2.7.3. Small-world parameter

To estimate small-worldness of the developing networks, we applied the formalism by Humphries and Gurney (2008):

$\begin{array}{cc}s=\frac{C/C^{rand}}{L/L^{rand}}& \left(11\right)\end{array}$

We replaced the clustering coefficient C and the characteristic path length L by the corresponding versions for weighted directed graphs as described above. C^rand and L^rand were taken from an Erdős-Rényi random graph generated with the same number of synapses as in the developing networks at every update in connectivity.

2.7.4. Betweenness centrality

The local betweenness centrality is a measure for the importance of a node in a network. A high betweenness centrality of a node i means that many shortest paths between any two nodes k and l pass through node i. Thus, the local betweenness centrality counts the number of times, σ_kl(i), that node i is on a shortest path between two nodes k and l. The local betweenness centrality is normalized by the number σ_kl of alternative shortest paths between k and l that do not pass through node i. Global betweenness centrality is the sum of all local betweenness centrality values of the individual nodes:

$\begin{array}{cc}B{C}_{global}={\displaystyle \sum _i^{n^{ex}}{\displaystyle \sum _{k\ne i\ne l}\frac{{\sigma }_{kl}(i)}{{\sigma }_{kl}}}}& \left(12\right)\end{array}$

Consequently, betweenness centrality provides a measure for how well networks are interconnected. A high local or global betweenness centrality means that individual nodes or the entire network, respectively, is badly interconnected, because all information has to travel through the same nodes in the absence of alternative routes or by-passes. Clinical data shows that after brain lesions, betweenness centrality of directly and indirectly affected brain areas changes (Wang et al., 2010). Note that as for the measurements above, the shortest paths are based on weighted excitatory connections W^ex,ex_i,j. Therefore, global betweenness centrality was computed by the formalism for weighted directed networks by Brandes (2001) as implemented in the brain connectivity toolbox (Rubinov and Sporns, 2010).

2.7.5. Efficiency

Global efficiency is related to the inverse characteristic path length with the advantage over characteristic path length that it can be meaningfully computed also of disconnected graphs (Latora and Marchiori, 2001; Achard and Bullmore, 2007). While path lengths between disconnected cells are infinite, efficiency becomes zero and, thus, adds neutrally to global efficiency.

$\begin{array}{cc}{E}_{global}=\frac{1}{n^{ex}\left(n^{ex}-1\right)}{\displaystyle \sum _{i\ne j}\frac{1}{d_{i,j}}}& \left(13\right)\end{array}$

where n^ex is the number of excitatory neurons.

2.7.6. Euclidean distance

Although not a topology measure in a strict graph theoretical sense, the average Euclidean distance between nodes that are connected by synapses was measured in order to help interpret the development of characteristic path length and clustering coefficient. To obtain the average Euclidean distance ED, we multiplied the Euclidean distances between all pairs of excitatory neurons in the network with the number of synapses between these neurons. We then summed up all of the so-weighted distances and divided the sum by the number of excitatoy synapses:

$\begin{array}{cc}ED=\left({\displaystyle \sum _{i,j}^{n^{ex}}\sqrt{{\left(po{s}_{xj}-po{s}_{xi}\right)}^2+{\left(po{s}_{yj}-po{s}_{yi}\right)}^2}}{W}_{i,j}^{ex,ex}\right)/{\displaystyle \sum _{i,j}^{n^{ex}}{W}_{i,j}^{ex,ex}}& \left(14\right)\end{array}$

with i,j ∈ {Ex}

3. Results

We started each network simulation with zero connectivity and zero synaptic elements. Due to the homeostatic formation of axonal and dendritic elements, model neurons are able to form synapses and to adapt the number of synapses so as to reach a homeostatic equilibrium in electrical activity (Figure 2). For a wide range of set-points ϵ, model neurons adapt their average firing rate so that, at the end of each simulation (here after 15,000 updates in connectivity), their calcium concentration [Ca²⁺] has converged to ϵ, with a corresponding firing rate y (in Hz) that follows the linear relation y = 100 × [Ca²⁺] (Figure 3). We investigated how the topology of these self-organizing networks developed over time when neurons favor short-range connections over long-range connections, or, alternatively, when all connections are equally likely. In the first case, as we will show in detail later on, networks developed a distinct small-world property, with small-world parameter s markedly greater than 1, whereas in the second case, s reaches 1. Therefore, we will call networks that resulted from distance-dependent synapse formation small-world networks and those that resulted from synapse formation without distance-dependency random networks.

FIGURE 2

Figure 2. Development of intracellular calcium concentration [Ca²⁺] over time. (A) Development of calcium concentration in small-world networks arising from synapse formation that favors short-range over long-range connections. Mean calcium concentrations averaged over five simulations reach the set-point ϵ = 0.7 when the number of synaptic elements is regulated homeostatically (yellow), whereas calcium concentrations remain much lower when there is no homeostatic regulation (gray). (B) Development of calcium concentration in random networks without any distance-dependency in synapse formation. The random network with homeostatic regulation of synaptic elements (blue) also reaches the homeostatic set-point ϵ, whereas the network without homeostasis (gray) and the same number of synapses as the homeostatic network has much lower values in calcium. Shadow around the curves (hardly visible since so small) indicates the standard deviation.

FIGURE 3

Figure 3. Neurons develop their connectivity in order to reach a homeostatic set-point ϵ of intracellular calcium. The figure shows the firing rates attained for different set-point values of calcium. For each calcium concentration, the firing rates of all neurons were pooled from four different simulations recorded over the last 20,000 ms. The central mark of each box is the median firing rate; boxes represent the 25th and 75th percentiles; the whiskers extend to the most extreme data points not considered “outliers”; and “outliers” are plotted individually (which show up here as thick bold lines). An “outlier” is a value that is more than 1.5 times the interquartile range away from the top or bottom of the box. The firing rate is proportional to ϵ by a factor of 100. For this set of simulations lower external inputs I^ext with mean 2 mVms⁻¹ were used.

In networks that favor long-range over short-range connections (small-world networks) (Figure 4A), s constantly increased and reached a maximum markedly greater than 10 at around 7000 updates in connectivity. Small-world networks that were set-up by an additional rule for homeostasis in electrical activity reached a plateau of about s = 10 very early but began to decrease again around the time that neuronal activities reached the homeostatic set-point ϵ. Nevertheless, small-world networks with homeostasis maintained their small-world property throughout the whole course of the simulation (s > 5 at T = 15,000). In networks without distance-dependency in synapse formation (random networks) (Figure 4B), s equaled 1 from the very beginning of network development. Random networks with and without homeostasis did not differ in the course of s (Figure 4B). Hence, homeostasis did not seem to have an impact on the development of topology in random networks. Knowing that homeostasis influenced the development of small-worldness, we further analyzed how homeostasis exerted its influence during network formation. For some simulations, s could not be computed for the first few updates in connectivity because of division by zero.

FIGURE 4

Figure 4. Small-worldness of developing networks. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations.

3.1. Homeostasis Influences the Clustering Coefficient and Characteristic Path Length of Developing Small-World Networks

The mean or characteristic path length in small-world networks without homeostasis showed, after an initial sharp rise, a steady decrease and converged toward values of around 3 (Figure 5A). Homeostatic small-world networks also started with a sharp rise and a subsequent decrease in characteristic path length. The decrease was even steeper than in the non-homeostatic case, and the characteristic path length converged toward slightly lower values than in non-homeostatic networks. In addition, networks with homeostasis showed a second but minor decrease in characteristic path length when activities reached the homeostatic set-point ϵ. The final values of the characteristic path length after 15,000 updates in connectivity were stable in both homeostatic and non-homeostatic small-world networks. The initial sharp increase is caused by the limited number of synapses in the network at the beginning of the simulation. Characteristic path lengths in developing random networks showed an identical course with and without homeostasis (Figure 5B). The final values in random networks were marginally lower than those in homeostatic small-world networks.

FIGURE 5

Figure 5. Characteristic path length in developing networks. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. The inset in (A) is a close-up of the time interval from 5000 to 10,000 updates in connectivity that clearly shows the decay in characteristic path length in homeostatic networks compared to non-homeostatic networks. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations.

The clustering coefficient in developing homeostatic small-world networks took a considerably different course from the coefficient in small-world networks without homeostasis (Figure 6A). Whereas the clustering coefficient in networks without homeostasis converged, with a small overshoot, toward high levels of over 1.6, networks with homeostasis generated a considerable overshoot and ended up at much lower values as compared with non-homeostatic networks. After a ramp-up phase, homeostatic networks reached a maximum clustering coefficient of about one; thereafter, the clustering coefficient decreased again with a decreasing negative slope. The maximum clustering coefficient was reached when the average electrical activities approached the homeostatic set-point ϵ. By contrast, we did not see an overshoot in clustering coefficient in homeostatic and non-homeostatic random networks (Figure 6B). Therefore, the homeostasis in electrical activity had no effect on clustering in random networks.

FIGURE 6

Figure 6. Clustering coefficient in developing networks. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations.

3.2. Homeostasis Favors Long-Range Connections in Small-World Networks

Particularly the development of the clustering coefficient, with its pronounced overshoot, determines the emerging small-worldness in networks favoring short-range connections over long-range connections. We therefore further studied how homeostasis influenced clustering in small-world networks. We tested the hypothesis that homeostasis produced more long-range connections than expected from the kernel function K (Equation 5). Computing the average Euclidean distance between the pre- and postsynaptic neuron for every synapse indeed revealed longer average Euclidean distances in homeostatic than in non-homeostatic networks (Figure 7A). In non-homeostatic networks, the average Euclidean distance was constant, because it is directly derived from the kernel function (in Equation 5: σ = 150 μm). In homeostatic networks, however, we observed two different phases. First, during initial network development ([Ca²⁺] << ϵ), the average Euclidean distance converged quickly toward a stable value of around 2, which was only slightly higher than in non-homeostatic networks. Secondly, when calcium concentrations approached the homeostatic set-point ϵ, the average Euclidean distances ramped up and reached values greater than 4. The considerable increase in the average Euclidean distance of synaptic connections coincided with a drop in clustering coefficient. As the initial high clustering of the network is due to the kernel function for synapse formation favoring short- over long-range connections, we may conclude that increasing Euclidean lengths of synaptic connections give rise to the decrease in clustering right at the time neurons approach the homeostatic set-point in electrical activity. The effect was also noticeable in the course of the characteristic path length but much less pronounced. We did not observe a comparable effect of homeostasis on Euclidean distances of synaptic connections in random networks (Figure 7B).

FIGURE 7

Figure 7. Average Euclidean distance between pre- and postsynaptic neurons for every synapse during network development. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations. Inset in (A) depicts the change in Euclidean distances for homeostatic random (blue) and homeostatic small-world networks (orange) with torus boundary conditions.

Are the increasing average Euclidean distances of synaptic connections caused by the fact that average neuronal activities are reaching a homeostatic equilibrium, or is this just some network effect that merely coincides with neurons equilibrating their activities? To answer this question, we assessed the number of all types of vacant elements (i.e., excitatory or inhibitory axonal elements, excitatory or inhibitory dendritic elements) and the spatial position of their hosting neurons at every time step of the simulation. We first checked whether there was any bias in the spatial position of synaptic elements that could generate more distant connections, e.g., a placement of synaptic elements at the boundaries of the network. We accumulated the number of vacant axonal (Figure 8A) and dendritic elements (Figure 8B) for each neuron for the first 5000 updates in connectivity with [Ca²⁺] < ϵ as well as for the next 5000 time steps with [Ca²⁺] ≈ ϵ. In the beginning of the simulation, we indeed found a little more vacant synaptic elements at the network boundaries, which can be explained by the fact that neurons at boundaries have less neighbors than neurons in the center and compensate for this by getting a higher number of vacant dendritic elements. However, by the time that the homeostatic equilibrium is reached, vacant synaptic elements were equally distributed over the network, and therefore the placement of vacant synaptic elements cannot account for the increasing distances of synaptic connections once the network has reached the homeostatic set-point. We found a comparable course of Euclidean distances in homeostatic small-world networks with torus boundary conditions (see inset of Figure 7A).

FIGURE 8

Figure 8. Spatial position of vacant synaptic elements in small-world networks. The purpose of this figure is to rule out that the position of vacant synaptic elements alone gave rise to longer-ranged synaptic connections. The left column shows vacant synaptic elements accumulated over the first 5000 updates in connectivity, whereas the right column shows the accumulation of vacant synaptic elements over the following 5000 updates. Although some tendency of vacant synaptic elements being located at the borders of the network is visible, this effect is gone after 5000 updates in connectivity. So there is no bias in the distribution of vacant synaptic elements when activity reaches the homeostatic set-point, and therefore the position of vacant synaptic elements alone cannot account for the formation of longer-ranged synaptic connections for T > 5000. (A) Vacant axonal elements. (B) Vacant dendritic elements. Color scale indicates number of vacant synaptic elements per neuron.

Since we could exclude a spatial bias in the distribution of vacant synaptic elements, we tested whether the increasing distances of synaptic connections were a direct consequence of the activity-dependent growth rules. On the basis of the number of vacant elements per neuron and the Euclidean distance between the host neurons, we determined the most likely synaptic connection every time a new synapses was formed, which is equivalent to the maximum of P^form_i,j (Equation 7). Since the change in the number of synaptic elements and therefore also the distribution of vacant synaptic elements are activity-dependent, the most likely synaptic connection to be formed is a direct consequence of the neurons' electrical activities. It turned out that in the beginning of development, when all neurons offer vacant synaptic elements, the most likely synaptic connections are those between adjacent neurons (Figures 9A,B). In other words, when neuronal activities were much lower than the homeostatic set-point ϵ, the kernel function had a large impact on the Euclidean length of synaptic connections. Therefore, at this stage of development, homeostatic and non-homeostatic networks did not differ much in connection lengths. However, when the activity of all individual neurons approached the set-point ϵ, vacant synaptic elements became rare and matching synaptic elements were available, if at all, only between more distant neurons. Expected distances for most likely synaptic connections therefore became much larger. Moreover, the distribution of expected Euclidean lengths of new synaptic connections did not follow the Gaussian-shaped kernel any longer but became much wider and flatter. At the same time, it took much longer for a synapse to form because, although longer-range synaptic connections were the most likely ones, their probability was still very low due to the kernel function. Due to a lack of shorter-range alternatives, these longer-range synaptic connections were nonetheless formed at some point in time, because the kernel function is non-zero for all distances. Taken together, the Euclidean length of synaptic connections in small-world networks is influenced by the homeostatic formation of synaptic elements. As expected from the topology data, no activity-dependent effect on the Euclidean distance of synaptic connections was observed in random networks.

FIGURE 9

Figure 9. The spatial distribution of newly formed synapses change in dependence on the calcium concentration. In (A,B), each dot represents the Euclidean distances between those neurons that are most likely to form a synaptic connection with each other at this update in connectivity. For this, we took at every update in connectivity in which vacant synaptic elements were available the Euclidean distance of the connection from neuron j to i for which P^form_i,j (Equation 7) was maximal. In (C,D), each dot represents the length of that connection (again in terms of the Euclidean distance between the connected neurons) for which synapse deletion was most likely, i.e., P^del_i,j (Equation 6) was maximal for every update in connectivity in which synapses had to be deleted. In (A,C), we plotted the Euclidean distances for synapse formation and deletion over time. The black curve (right y-axis) indicates the course of the calcium concentration [Ca²⁺]. In (B,D), we plotted synapse formation and deletion in dependence on [Ca²⁺]. The color code in all panels indicates the density of the dots in the diagrams, with blue and red representing low and high densities of dots, respectively. The figure essentially shows that before calcium reaches the homeostatic set-point ϵ, the distribution for synapse formation is rather Gaussian, following the Kernel function K (Equation 5). The distribution becomes random and scattered, with increased Euclidean distances, when calcium is at the set-point. The stripes in the distribution arise from the fact that not all Euclidean distances are possible due to the grid layout of the network. There is no change in the distribution for synapse deletion.

Additionally, synapse deletion could in principle also influence the Euclidean length of synaptic connections if, for example, with increasing neuronal activities preferentially short connections would be deleted. Therefore, we further tested whether the expected Euclidean length of synaptic connections that were most likely to be deleted correlated with the current average activity in the network. However, an activity-dependent effect on synapse deletion was not observed (Figures 9C,D).

3.3. Homeostasis Decreases the Betweenness Centrality in Small-World Networks

Betweenness centrality (or inbetweenness) is a measure for the importance of nodes in a network. Compared with random networks, small-world networks without homeostatic synapse formation had a relatively high betweenness centrality (Figure 10A). Small-world networks with homeostasis, by contrast, revealed a pronounced decrease in betweenness centrality over developmental time. Values peaked in the beginning of development when neurons first connected to their nearest neighbors only. However, right after the peak, homeostatic synapse formation generated topologies with values for betweenness centrality that were lower than in non-homeostatic networks. Betweenness centrality is high in non-homeostatic small-world networks because synapse formation is only determined by the kernel function (Equation 5), which often caused that the same cell pairs were connected repeatedly. The consequence is a limited number of shortest paths between any pairs of nodes in the network. By contrast, in homeostatic networks, synapse formation depends on the availability of synaptic elements, which can force synapses to be formed that are less likely according to the kernel. This may increase the number of alternative paths between two neurons and therefore also increase the number of multiple shortest paths through the network. Hence, betweenness centrality quickly decreased over time.

FIGURE 10

Figure 10. Changing betweenness centrality over time. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations.

Because the variety of shortest paths in non-homeostatic networks is limited, betweenness centrality reached a stable plateau in these networks. Interestingly, the betweenness centrality in homeostatic networks initially also converged toward a quasi-stable level. However, as soon as activities approached the homeostatic set-point ϵ, betweenness centrality strongly decreased. Over time, the rate of decrease slowed down. The decrease in betweenness centrality precisely coincided with the increase in Euclidean lengths of synaptic connections. Because the kernel function favors short- over long-range connections and therefore initially creates networks with high betweenness centrality, we may conclude that the formation of longer-range connections (in an Euclidean sense) lead to a decreasing betweenness centrality. Any new long-range connection not present in the network before creates new shortest paths, which in turn decreases betweenness centrality. By contrast, homeostasis did not affect the course of betweenness centrality in random networks (Figure 10B).

3.4. Small-World Networks Become more Efficient by Homeostasis

Efficient information transmission is probably the most-desired property in computational networks. Small-world networks are very efficient because they combine a high clustering coefficient with a short characteristic path length. Nevertheless, their efficiency is still markedly lower than that of random networks. Homeostatic small-world networks, however, generated efficiency levels that during the whole course of development exceeded the levels in non-homeostatic small-world networks (Figure 11A). Remarkably, at the time when average activities reached the homeostatic set-point ϵ, efficiency levels of homeostatic networks further increased and almost reached the levels in random networks (Figure 11B). Consequently, favoring more distant connections in combination with homeostasis in electrical activity led to a more efficient network topology than achieved without homeostasis.

FIGURE 11

Figure 11. Global efficiency changes over time due to network development. (A) Small-world networks with (orange) and without (gray) homeostasis in electrical activity. (B) Random networks with (blue) and without (gray) homeostasis in electrical activity. Means over five simulations per scenario. Shadings of the curves indicate standard deviations.

4. Discussion

We have shown that network formation favoring short-range over long-range connections produced networks with a pronounced small-world structure. Networks with homeostasis in electrical activity developed a weaker small-world structure in favor of more efficient wiring of connections. Global efficiency particularly increased when network activity reached the homeostatic set-point. Increased global efficiency was caused by the fact that homeostasis favored longer-ranged connections, which affected clustering as well as characteristic path length. Thus, network topology continued to change even after the network had reached a homeostatic equilibrium in electrical activity.

Adding more long-range connections to a small-world network makes the network more efficient but also more random. This is apparent in our simulations, too, by a decreasing clustering coefficient and a decreasing characteristic path length. Nevertheless, the small-world property of the networks was preserved throughout the whole course of development, and the decrease in clustering coefficient was slowing down over time. However, networks would most likely turn into random networks if rewiring continued indefinitely. Consequently, there seems to be a trade-off in network development between high clustering and strong small-worldness on the one hand and more randomness and higher efficiency on the other hand. The latter particularly arises when networks continue to rewire their circuitry when they are in a homeostatic equilibrium.

With additional long-range connections betweenness centrality decreases. Networks with low betweenness centrality are more robust against lesions because all nodes are equally important. By contrast, networks with high betweenness centrality are very vulnerable to lesions. If neurons that are part of many shortest paths are lost, the characteristic path length will immediately increase, with a significant impact on information transmission. It is remarkable that a self-organizing process that forms networks by striving toward homeostasis in electrical activity as a side effect produces topologies with lower betweenness centrality that contribute to a higher robustness against lesions.

Key to the homeostasis-driven change in topology is the increasing Euclidean length of connections. Since we did not observe an increase in connection length in networks without homeostasis, we concluded that the increase in connection lengths was caused by homeostasis in electrical activity. To create networks without homeostasis, we took the number of synapses from the homeostatic network at every update in connectivity and distributed them randomly under the same kernel in the non-homeostatic network. There are, in fact, other ways to create networks without homeostasis. We additionally built non-homeostatic networks by initially giving all model neurons a fixed number of vacant synaptic elements and then running updates in connectivity until no more synapses could be formed. Also in this scenario we did not see an increase in connection length. In a third scenario, we added a few vacant synaptic elements to all neurons before every update in connectivity. No matter how few vacant elements we added and how long we ran the network, we did not obtain more long-range connections than expected from the kernel. We only observed more long-range connections when we slowly added vacant elements at a few randomly selected neurons at a time (i.e., spatial sparseness in the formation of synaptic elements). From these observations we concluded that the contribution of homeostasis is not only to limit the number of available vacant synaptic elements but also to generate a certain sparseness in the formation of vacant elements. In fact, as long as the neurons were far away from the homeostatic set-point, synaptic element formation was not sparse at all as all neurons added elements roughly at the same time at equal rates. Only when the neurons approached the homeostatic set-point did sparseness in synaptic element formation arise. Homeostasis creates this sparseness because balancing the activity of one neuron immediately affects the activity in other neurons, which in turn may be driven further away from the homeostatic set-point and then start forming new vacant elements. The presence of inhibitory neurons would further reinforce this process.

Homeostasis in electrical activity is one way by which a local process can give rise to a change in global topology. Another comparable mechanisms was provided by Kaiser et al. (2009). In their model study, they showed that a simple axonal growth process can generate the experimentally observed exponential decrease in number of connections with increasing connection length. As a result of the growth process, most connections become short-range, although long-range connections also arise, but in lower numbers. The idea of their model is that axons grow out until they hit a postsynaptic target. The capacity of model neurons to receive connections is limited, and if nearby target neurons are completely occupied by incoming connections, axons continue to grow out until they hit a vacant target. Hence, Euclidean connection length increases over developmental time. In this model, the spatial growth process in combination with a hard boundary on the number of incoming synapses per model neuron generates the increase in connection lengths.

There are striking similarities between the topology of our model networks and that of developing dissociated cell cultures. It is well known that cultured neuronal networks can form small-world topology (Bettencourt et al., 2007; Yu et al., 2008; Gerhard et al., 2011). Downes et al. (2012) reported an increase in clustering coefficient of dissociated cell cultures between 14 days in vitro (DIV) and 28 DIV and a subsequent drop until 35 DIV, a course of development that is comparable to the course of development in our model networks with homeostasis. Between 14 and 35 DIV, the mean shortest path length did not change significantly, only showing a slight drop around 28 DIV. Consequently, small-worldness reached its maximum around 28 DIV. Moreover, the experimental data indicated the presence of longer synaptic connections from 28 DIV onwards that had not been not present at 21 DIV. From our previous studies we know that dissociated cell cultures reach homeostasis around this time (Tetzlaff et al., 2010). Therefore, we may hypothesize that the increase in connection length in dissociated cell cultures may be due to neurons reaching a homeostatic equilibrium in electrical activity. Other synaptic plasticity mechanisms not currently incorporated in our model may of course also have contributed to network formation in developing cell cultures.

On a macroscopic scale, functional imaging data reveal a development of small-world topology that also has interesting similarities with the self-organizing network formation in our model. The infant human brain has small-world properties already at the third post-natal week (Fransson et al., 2011). During the following 2 years, network topology undergoes a significant refinement: brain networks increase their small-worldness, global efficiency and number of long-distance connections (Gao et al., 2011). Although our network model shows only an initial transient increase in small-worldness, it may offer a simple explanation for the sudden increase in number of long-range connections and the associated increase in global efficiency. Could it be that even in the human brain, neurons establishing a homeostatic equilibrium in electrical activity produce—as an emergent property of the homeostatic growth process—more long-range connections? Remarkably, the increase in number of long-range connections occurs not during a genetically-encoded formation of an initial embryonic layout of projections but during the post-natal critical period (Gao et al., 2011), during which neurons are highly sensitive to afferent input. In general, the importance of critical periods is to balance excitatory and inhibitory circuits and to establish homeostasis in neuronal electrical activity (Hensch, 2005; Butz et al., 2009b). Considering that long-range connections arise in local as well as global networks, our study raises the interesting hypothesis that homeostasis in electrical activity may be the driving force for the formation of long-range connections on both a microscopic and a macroscopic scale.

Homeostasis in electrical activity is a ubiquitous principle in the nervous system (Wolff and Wagner, 1983; Ramakers et al., 1991; Abbott and Nelson, 2000) and a variety of plasticity mechanisms can act homeostatically. Scaling of synaptic strengths (Turrigiano, 1999), for example, has been reported as a mechanism acting at existing synapses to stabilize postsynaptic firing in cortical, hippocampal and spinal cord networks (Lissin et al., 1998; O'Brien et al., 1998; Turrigiano and Nelson, 1998). Even in the mature brain, not only the strength of synapses but also the formation of new synapses can contribute to the stabilization of neuronal activity, for example after focal retinal lesions (Butz and van Ooyen, 2013). In developing dissociated cell cultures, we showed that homeostasis in electrical activity may be a precondition for the emergence of self-organized criticality in neuronal firing (Tetzlaff et al., 2010). In another study, we showed that homeostasis in electrical activity can regulate the synaptic embedding of newly formed neurons in the mature hippocampal dentate gyrus (adult neurogenesis) and can account for the experimentally observed counter-intuitive inverse relationship between cell proliferation rate and synaptogenesis (Butz et al., 2008).

In summary, we conclude that homeostatic regulation of electrical activity together with simple distance-dependent formation of connections is capable of creating, in a self-organizing manner, neuronal networks that are more robust and more efficient than networks grown without homeostatic regulation. Strikingly, the growth process revealed features of developing topologies that are also observed in dissociated cell cultures and infant human brains. The formation of network topology by a self-organizing, local growth process may also be relevant for automatically generating the connectivity structure of large-scale neuronal networks (Potjans and Diesmann, 2014) that are currently studied in enterprizes such as the Human Brain Project (www.humanbrainproject.eu).

Funding

This work was partly funded by the Helmholtz Association through the Helmholtz Portfolio Theme “Supercomputing and Modeling for the Human Brain.” Markus Butz and Ines D. Steenbuck were also partly supported by the NETFORM project (grant number 635.100.017, awarded to Arjen van Ooyen) of the Computational Life Sciences program of the Netherlands Organization for Scientific Research (NWO) and by a starting grant awarded to Markus Butz from the collaborative research center (SFB 874) funded by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG).

Conflict of Interest Statement

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

References

Abbott, L. F., and Nelson, S. B. (2000). Synaptic plasticity: taming the beast. Nat. Neurosci. 3, 1178–1183. doi: 10.1038/78829

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Abbott, L. F., and Rohrkemper, R. (2007). A simple growth model constructs critical avalanche networks. Prog. Brain Res. 165, 13–19. doi: 10.1016/S0079-6123(06)65002-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Achard, S., and Bullmore, E. (2007). Efficiency and cost of economical brain functional networks. PLoS Comput. Biol. 3:e17. doi: 10.1371/journal.pcbi.0030017

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bassett, D. S., and Bullmore, E. (2006). Small-world brain networks. Neuroscientist 12, 512–523. doi: 10.1177/1073858406293182

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bettencourt, L. M., Stephens, G. J., Ham, M. I., and Gross, G. W. (2007). Functional structure of cortical neuronal networks grown in vitro. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75, 021915. doi: 10.1103/PhysRevE.75.021915

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Binzegger, T., Douglas, R. J., and Martin, K. A. (2004). A quantitative map of the circuit of cat primary visual cortex. J. Neurosci. 24, 8441–8453. doi: 10.1523/JNEUROSCI.1400-04.2004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Borisyuk, R., Cooke, T., and Roberts, A. (2008). Stochasticity and functionality of neural systems: mathematical modelling of axon growth in the spinal cord of tadpole. Biosystems 93, 101–114. doi: 10.1016/j.biosystems.2008.03.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Braitenberg, V., and Schüz, A. (1998). Cortex: Statistics and Geometry of Neuronal Connectivity, 2nd edn. Berlin: Springer. doi: 10.1007/978-3-662-03733-1

CrossRef Full Text

Brandes, U. (2001). A Faster algorithm for betweenness centrality. J. Math. Soc. 25, 163–177. doi: 10.1080/0022250X.2001.9990249

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Butz, M., and Teuchert-Noodt, G. (2006). A simulation model for compensatory plasticity in the prefrontal cortex inducing a cortico-cortical dysconnection in early brain development. J. Neural Transm. 113, 695–710. doi: 10.1007/s00702-005-0403-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Butz, M., Wörgötter, F., and van Ooyen, A. (2009b). Activity-dependent structural plasticity. Brain Res. Rev. 60, 287–305. doi: 10.1016/j.brainresrev.2008.12.023

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cromme, L. J., and Dammasch, I. E. (1989). Compensation type algorithms for neural nets: stability and convergence. J. Math. Biol. 27, 327–340. doi: 10.1007/BF00275816

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dalva, M. B., Ghosh, A., and Shatz, C. J. (1994). Independent control of dendritic and axonal form in the developing lateral geniculate nucleus. J. Neurosci. 14, 3588–3602.

Pubmed Abstract | Pubmed Full Text

Dammasch, I. E., Wagner, G. P., and Wolff, J. R. (1986). Self-stabilization of neuronal networks. I. The compensation algorithm for synaptogenesis. Biol. Cybern. 54, 211–222. doi: 10.1007/BF00318417

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dammasch, I. E., Wagner, G. P., and Wolff, J. R. (1988). Self-stabilization of neuronal networks. II. Stability conditions for synaptogenesis. Biol. Cybern. 58, 149–158. doi: 10.1007/BF00364134

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Downes, J. H., Hammond, M. W., Xydas, D., Spencer, M. C., Becerra, V. M., Warwick, K., et al. (2012). Emergence of a small-world functional network in cultured neurons. PLoS Comput. Biol. 8:e1002522. doi: 10.1371/journal.pcbi.1002522

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fagiolo, G. (2007). Clustering in complex directed networks. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76:026107. doi: 10.1103/PhysRevE.76.026107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fransson, P., Aden, U., Blennow, M., and Lagercrantz, H. (2011). The functional architecture of the infant brain as revealed by resting-state fMRI. Cereb. Cortex 21, 145–154. doi: 10.1093/cercor/bhq071

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Franze, K. (2013). The mechanical control of nervous system development. Development 140, 3069–3077. doi: 10.1242/dev.079145

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gao, W., Gilmore, J. H., Giovanello, K. S., Smith, J. K., Shen, D., Zhu, H., et al. (2011). Temporal and spatial evolution of brain network topology during the first two years of life. PLoS ONE 6:e25278. doi: 10.1371/journal.pone.0025278

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gerhard, F., Pipa, G., Lima, B., Neuenschwander, S., and Gerstner, W. (2011). Extraction of network topology from multi-electrode recordings: is there a small-world effect? Front. Comput. Neurosci. 5:4. doi: 10.3389/fncom.2011.00004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Groc, L., Petanjek, Z., Gustafsson, B., Ben-Ari, Y., Hanse, E., and Khazipov, R. (2002). In vivo blockade of neural activity alters dendritic development of neonatal CA1 pyramidal cells. Eur. J. Neurosci. 16, 1931–1938. doi: 10.1046/j.1460-9568.2002.02264.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hensch, T. K. (2005). Critical period plasticity in local cortical circuits. Nat. Rev. Neurosci. 6, 877–888. doi: 10.1038/nrn1787

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hilgetag, C. C., and Kaiser, M. (2004). Clustered organization of cortical connectivity. Neuroinformatics 2, 353–360. doi: 10.1385/NI:2:3:353

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Humphries, M. D., and Gurney, K. (2008). Network ‘small-world-ness’: a quantitative method for determining canonical network equivalence. PLoS ONE 3:e0002051. doi: 10.1371/journal.pone.0002051

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hutchins, B. I., and Kalil, K. (2008). Differential outgrowth of axons and their branches is regulated by localized calcium transients. J. Neurosci. 28, 143–153. doi: 10.1523/JNEUROSCI.4548-07.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jourdain, P., Fukunaga, K., and Muller, D. (2003). Calcium/calmodulin-dependent protein kinase II contributes to activity-dependent filopodia growth and spine formation. J. Neurosci. 23, 10645–10649.

Pubmed Abstract | Pubmed Full Text

Kaiser, M., and Hilgetag, C. C. (2006). Nonoptimal component placement, but short processing paths, due to long-distance projections in neural systems. PLoS Comput. Biol. 2:e95. doi: 10.1371/journal.pcbi.0020095

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kaiser, M., Hilgetag, C. C., and van Ooyen, A. (2009). A simple rule for axon outgrowth and synaptic competition generates realistic connection lengths and filling fractions. Cereb. Cortex 19, 3001–3010. doi: 10.1093/cercor/bhp071

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Keck, T., Keller, G. B., Jacobsen, R. I., Eysel, U. T., Bonhoeffer, T., and Hubener, M. (2013). Synaptic scaling and homeostatic plasticity in the mouse visual cortex in vivo. Neuron 80, 327–334. doi: 10.1016/j.neuron.2013.08.018

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Keck, T., Mrsic-Flogel, T. D., Vaz Afonso, M., Eysel, U. T., Bonhoeffer, T., and Hubener, M. (2008). Massive restructuring of neuronal circuits during functional reorganization of adult visual cortex. Nat. Neurosci. 11, 1162–1167. doi: 10.1038/nn.2181

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kirov, S. A., Goddard, C. A., and Harris, K. M. (2004). Age-dependence in the homeostatic upregulation of hippocampal dendritic spine number during blocked synaptic transmission. Neuropharmacology 47, 640–648. doi: 10.1016/j.neuropharm.2004.07.039

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ko, H., Cossell, L., Baragli, C., Antolik, J., Clopath, C., Hofer, S. B., et al. (2013). The emergence of functional microcircuits in visual cortex. Nature 496, 96–100. doi: 10.1038/nature12015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Konur, S., and Ghosh, A. (2005). Calcium signaling and the control of dendritic development. Neuron 46, 401–405. doi: 10.1016/j.neuron.2005.04.022

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Latora, V., and Marchiori, M. (2001). Efficient behavior of small-world networks. Phys. Rev. Lett. 87:198701. doi: 10.1103/PhysRevLett.87.198701

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lissin, D. V., Gomperts, S. N., Carroll, R. C., Christine, C. W., Kalman, D., Kitamura, M., et al. (1998). Activity differentially regulates the surface expression of synaptic AMPA and NMDA glutamate receptors. Proc. Natl. Acad. Sci. U.S.A. 95, 7097–7102. doi: 10.1073/pnas.95.12.7097

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lohmann, C., and Wong, R. O. (2005). Regulation of dendritic growth and plasticity by local and global calcium dynamics. Cell Calcium 37, 403–409. doi: 10.1016/j.ceca.2005.01.008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McAssey, M. P., Bijma, F., Tarigan, B., van Pelt, J., van Ooyen, A., and de Gunst, M. (2014). A morpho-density approach to estimating neural connectivity. PLoS ONE 9:e86526. doi: 10.1371/journal.pone.0086526

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McKinney, R. A., Capogna, M., Durr, R., Gahwiler, B. H., and Thompson, S. M. (1999a). Miniature synaptic events maintain dendritic spines via AMPA receptor activation. Nat. Neurosci. 2, 44–49. doi: 10.1038/4548

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McKinney, R. A., Luthi, A., Bandtlow, C. E., Gahwiler, B. H., and Thompson, S. M. (1999b). Selective glutamate receptor antagonists can induce or prevent axonal sprouting in rat hippocampal slice cultures. Proc. Natl. Acad. Sci. U.S.A. 96, 11631–11636. doi: 10.1073/pnas.96.20.11631

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Missler, M., Wolff, A., Merker, H. J., and Wolff, J. R. (1993). Pre- and postnatal development of the primary visual cortex of the common marmoset. II. Formation, remodelling, and elimination of synapses as overlapping processes. J. Comp. Neurol. 333, 53–67. doi: 10.1002/cne.903330105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

O'Brien, R. J., Kamboj, S., Ehlers, M. D., Rosen, K. R., Fischbach, G. D., and Huganir, R. L. (1998). Activity-dependent modulation of synaptic AMPA receptor accumulation. Neuron 21, 1067–1078. doi: 10.1016/S0896-6273(00)80624-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Park, H. J., and Friston, K. (2013). Structural and functional brain networks: from connections to cognition. Science 342, 1238411. doi: 10.1126/science.1238411

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Potjans, T. C., and Diesmann, M. (2014). The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cereb. Cortex 24, 785–806. doi: 10.1093/cercor/bhs358

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ramakers, G. J., Raadsheer, F. C., Corner, M. A., Ramaekers, F. C., and Van Leeuwen, F. W. (1991). Development of neurons and glial cells in cerebral cortex, cultured in the presence or absence of bioelectric activity: morphological observations. Eur. J. Neurosci. 3, 140–153. doi: 10.1111/j.1460-9568.1991.tb00074.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52, 1059–1069. doi: 10.1016/j.neuroimage.2009.10.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Scheiffele, P., Fan, J., Choih, J., Fetter, R., and Serafini, T. (2000). Neuroligin expressed in nonneuronal cells triggers presynaptic development in contacting axons. Cell 101, 657–669. doi: 10.1016/S0092-8674(00)80877-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Siegel, F., and Lohmann, C. (2013). Simultaneous imaging of structural plasticity and calcium dynamics in developing dendrites and axons. Cold Spring Harb. Protoc. 2013:pdb.prot078592. doi: 10.1101/pdb.prot078592

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O., and Zwi, J. D. (2004). The small world of the cerebral cortex. Neuroinformatics 2, 145–162. doi: 10.1385/NI:2:2:145

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stan, A., Pielarski, K. N., Brigadski, T., Wittenmayer, N., Fedorchenko, O., Gohla, A., et al. (2010). Essential cooperation of N-cadherin and neuroligin-1 in the transsynaptic control of vesicle accumulation. Proc. Natl. Acad. Sci. U.S.A. 107, 11116–11121. doi: 10.1073/pnas.0914233107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tailby, C., Wright, L. L., Metha, A. B., and Calford, M. B. (2005). Activity-dependent maintenance and growth of dendrites in adult cortex. Proc. Natl. Acad. Sci. U.S.A. 102, 4631–4636. doi: 10.1073/pnas.0402747102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tetzlaff, C., Okujeni, S., Egert, U., Wörgötter, F., and Butz, M. (2010). Self-organized criticality in developing neuronal networks. PLoS Comput. Biol. 6:e1001013. doi: 10.1371/journal.pcbi.1001013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turrigiano, G. G. (1999). Homeostatic plasticity in neuronal networks: the more things change, the more they stay the same. Trends Neurosci. 22, 221–227. doi: 10.1016/S0166-2236(98)01341-1

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turrigiano, G. G., and Nelson, S. B. (1998). Thinking globally, acting locally: AMPA receptor turnover and synaptic strength. Neuron 21, 933–935. doi: 10.1016/S0896-6273(00)80607-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Uesaka, N., Hirai, S., Maruyama, T., Ruthazer, E. S., and Yamamoto, N. (2005). Activity dependence of cortical axon branch formation: a morphological and electrophysiological study using organotypic slice cultures. J. Neurosci. 25, 1–9. doi: 10.1523/JNEUROSCI.3855-04.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Ooyen, A. (2011). Using theoretical models to analyse neural development. Nat. Rev. Neurosci. 12, 311–326. doi: 10.1038/nrn3031

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Ooyen, A., Carnell, A., de Ridder, S., Tarigan, B., Mansvelder, H. D., Bijma, F., et al. (2014). Independently outgrowing neurons and geometry-based synapse formation produce networks with realistic synaptic connectivity. PLoS ONE 9:e85858. doi: 10.1371/journal.pone.0085858

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Ooyen, A., and van Pelt, J. (1994). Activity-dependent neurite outgrowth and neural network development. Prog. Brain Res. 102:245–259. doi: 10.1016/S0079-6123(08)60544-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Oss, C., and van Ooyen, A. (1997). Effects of inhibition on neural network development through activity-dependent neurite outgrowth. J. Theor. Biol. 185, 263–280. doi: 10.1006/jtbi.1996.0361

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Pelt, J., and van Ooyen, A. (2013). Estimating neuronal connectivity from axonal and dendritic density fields. Front. Comput. Neurosci. 7:160. doi: 10.3389/fncom.2013.00160

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, L., Yu, C., Chen, H., Qin, W., He, Y., Fan, F., et al. (2010). Dynamic functional reorganization of the motor execution network after stroke. Brain 133, 1224–1238. doi: 10.1093/brain/awq043

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. Nature 393, 440–442. doi: 10.1038/30918

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wolff, J., and Wagner, G. (1983). “Selforganization in synaptogenesis: interaction between the formation of excitatory and inhibitory synapses,” in Synergetics of the Brain, eds E. Basar, H. Flohr, H. Haken, and A. Mandell (Berlin; Heidelberg; New York; Tokyo: Springer).

Wong, R. O., and Ghosh, A. (2002). Activity-dependent regulation of dendritic growth and patterning. Nat. Rev. Neurosci. 3, 803–812. doi: 10.1038/nrn941

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yamamoto, N., Tamada, A., and Murakami, F. (2002). Wiring of the brain by a range of guidance cues. Prog. Neurobiol. 68, 393–407. doi: 10.1016/S0301-0082(02)00129-6

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yu, S., Huang, D., Singer, W., and Nikolic, D. (2008). A small world of neuronal synchrony. Cereb. Cortex 18, 2891–2901. doi: 10.1093/cercor/bhn047

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text