Temporal ordering of input modulates connectivity formation in a developmental neuronal network model of the cortex

Preterm infant brain activity is discontinuous; bursts of activity recorded using EEG (electroencephalography), thought to be driven by subcortical regions, display scale free properties and exhibit a complex temporal ordering known as long-range temporal correlations (LRTCs). During brain development, activity-dependent mechanisms are essential for synaptic connectivity formation, and abolishing burst activity in animal models leads to weak disorganised synaptic connectivity. Moreover, synaptic pruning shares similar mechanisms to spike-timing dependent plasticity (STDP), suggesting that the timing of activity may play a critical role in connectivity formation. We investigated, in a computational model of leaky integrate-and-fire neurones, whether the temporal ordering of burst activity within an external driving input could modulate connectivity formation in the network. Connectivity evolved across the course of simulations using an approach analogous to STDP, from networks with initial random connectivity. Small-world connectivity and hub neurones emerged in the network structure—characteristic properties of mature brain networks. Notably, driving the network with an external input which exhibited LRTCs in the temporal ordering of burst activity facilitated the emergence of these network properties, increasing the speed with which they emerged compared with when the network was driven by the same input with the bursts randomly ordered in time. Moreover, the emergence of small-world properties was dependent on the strength of the LRTCs. These results suggest that the temporal ordering of burst activity could play an important role in synaptic connectivity formation and the emergence of small-world topology in the developing brain.


Introduction
Network connectivity shapes activity and modulates information transfer in the brain. For example, small-world network architecture allows e cient integration and segregation of information [1], and hub neurones or regions play a key role in carrying information throughout the brain [2]. This structured connectivity emerges during early human brain development; a small-world modular network organisation with hub nodes can be observed in preterm di↵usion and functional MRI, with a significant increase in small-world topology between 30 and 40 weeks' gestation [3,4].
The major period of connectivity formation and refinement in the cortex starts during foetal development from approximately 20 weeks' gestation and continues for the first few years of postnatal life [5,6]. MEG recordings of foetal brain activity and EEG recordings from preterm infants are characterised by discontinuous activity -bursts of slow wave oscillations with nested high frequency activity are interspersed within periods of apparent electrical silence [7].
These bursts can occur in response to sensory stimulation [8][9][10], or following movement [11], but the majority occur spontaneously in the background EEG [11]. Spontaneous bursts are thought to originate from regions such as the subplate [12,13], a transient population of neurones present in early development [14], and may also relate to activity in the insula [15].
Neuronal activity is crucial for connectivity formation [16], and blocking or reducing burst activity during critical developmental periods, for example, through removal of the subplate in animal models, leads to abnormal cortical network connectivity, with weak thalamocortical connectivity [17] and loss of cortical columnar structure [18,19].
The temporal organisation of neuronal activity in early development may also play a key role in connectivity formation -rearing fish in an environment with stroboscopic illumination prevents the refinement of retinotectal maps [20] and periodic electrical stimulation of the ferret optic nerve results in altered orientation selectivity in the cortex [21]. Moreover, activity-dependent mechanisms at the synaptic level are key to connectivity refinement in the developing brain, and synaptic pruning shares similar molecular pathways with long-term depression (LTD) in the adult brain [22], suggesting that the temporal organisation of activity may play a critical role in connectivity formation. Recently it has been demonstrated that the burst activity in preterm EEG exhibits scale-free properties [23], and that the bursts do not occur randomly in time but follow a complex temporal ordering, known as long-range temporal correlations (LRTCs) [24]. Thus, in the preterm neonatal brain, the timing of any given EEG burst is correlated with the time of occurrence of all previous bursts of EEG activity [24]. Whether this correlated temporal structure of the timing of burst activity a↵ects connectivity formation in the developing brain is an important open question.
Here we consider this question by investigating connectivity formation in a simple activitydependent neuronal network model of the cortex. Motivated by the excitatory role of GABA within the developing brain [25], we consider a model where all connections are excitatory.
We assume the cortex is driven by bursts of activity from a non-cortical source such as the subplate [12,13]. We compare connectivity formation in a network driven by bursts which exhibit LRTCs, to the network connectivity that emerges when the network is driven by bursts with random temporal ordering. This preserves the inter-burst interval distribution whilst changing the temporal correlations in the ordering of the inter-burst intervals. We also compare with the connectivity in networks driven by periodic bursts. Finally, we investigate the relationship between network connectivity parameters and the strength of LRTCs within the burst temporal organisation. We test the hypothesis that LRTCs in the burst activity of the external input promotes the emergence of small-world connectivity in the developing cortex.

Results
We investigated connectivity formation in directed networks of leaky integrate-and-fire neurones. Whilst this work is motivated by the observation of LRTCs in the EEG of preterm infants, i.e. an observation at the macroscopic scale, LRTCs have been demonstrated at multiple levels of the nervous system in adults, including in the spontaneous spiking of individual neurones [26]. Although the external input which provides the LRTCs within our model is on the macroscopic scale, with bursts originating from non-cortical sources such as the subplate [12,13], our focus here is on the e↵ect such macroscopic phenomenon has on connectivity formation at the microscopic (neuronal) scale.
Individuals neurones in our model received both this external input and input from connected neurones when these neurones fired. To ensure that neuronal firing did not saturate the network or die out completely, synaptic weights were updated according to the number of connections within the network (see Methods), which can be considered a form of homeostatic plasticity [27]. The external input to the network was a bursty input, Fig. 1A. In all cases the bursts themselves were of a fixed duration and amplitude; di↵erences in the driving input were only reflected in the temporal ordering of the bursts, i.e, in the ordering of the interburst intervals (IBIs) -the time between the bursts. In the case where the network was driven with burst dynamics that exhibit LRTCs, the sequence of IBIs exhibited long-range temporal correlations with a Hurst exponent (H) greater than 0.5. In the shu✏ed input case, the same IBIs were randomly re-ordered in time, giving a Hurst exponent of the sequence of IBIs of H ⇡ 0.5. Thus, in both cases the external input had the same IBI distribution, but the temporal ordering of the bursts was changed and in the latter case the temporal correlations were lost.
Driving the neuronal networks with bursty input also led to bursts within the network ( Fig. 1B,E). The Hurst exponent for the sequence of IBIs was determined using detrended fluctuation analysis (DFA). The DFA exponents for the sequences of IBIs in the network firing dynamics reflected those of the corresponding DFA exponents of the external input ( Fig. 1C,F), indicating that the temporal characteristics of the driving input were transmitted to the network activity itself.

Emergence of small-world topology and hub neurones
Network connectivity was allowed to evolve during simulations using an approach analogous to spike-timing dependent plasticity (STDP). The likelihood of gaining a connection between any two neurones was increased if the presynaptic neurone frequently fired just before the  connections (see Methods). Network connectivity was initially random, with 40% connectivity in a network with 200 neurones (alterations to these parameters are considered in the section Varying network size and density). As with STDP in the adult brain [28], and as similar mechanisms to LTD play a dominate role in neuronal network development [22], we initially set depression to be slightly stronger than potentiation. We also modelled alternatives, i.e. equal amounts of depression and potentiation and stronger potentiation than depression, which are described below (see section Varying the plasticity parameters).
With slightly stronger depression than potentiation, and driven with burst dynamics which exhibit LRTCs, the network on average lost connections across the course of the simulation and, although the e↵ect size is relatively small, the normalised clustering coe cient and small-world index increased. At the end of the simulations, the node in degree distribution is skewed, with some neurones showing much higher degree than others, indicating the presence of hub neurones (Fig. 2D).
Thus, across the course of the simulation a small-world topology and hub neurones emerge.
Moreover, the speed at which the proportion of connections, normalised clustering coe cient and small-world index changes is high at the start of the simulations, indicating that the network rapidly evolves to have connectivity with small-world properties.
In contrast, when the network is driven with the same external input but in which the burst order has been randomly shu✏ed in time (and therefore does not exhibit LRTCs but instead H ⇡ 0.5), the rate of change in the network parameters at the start of the simulation is much lower than with external input which exhibits LRTCs, Fig   The di↵erence in emergence of small-world properties when the external input exhibits LRTCs leads to the question as to whether the magnitude of the LRTCs also a↵ects connectivity formation. Therefore, we next investigated how driving the network with burst activity with IBIs which exhibit di↵erent strengths of LRTCs alters the rate of emergence of small-world properties.

Varying the Hurst exponent modulates connectivity formation
We investigated driving the network with burst dynamics which had 4 di↵erent levels of Hurst exponent: H ⇡ 0.5, 0.6, 0.7 and 0.8. In all cases, on average the network lost connections, but networks driven by external input with higher Hurst exponents lost connections at a faster initial rate and quickly evolved to a network with small-world properties, see Fig. 3. The average Hurst exponent estimated using DFA in preterm infants is 0.68, with a range of 0.55 -0.81 [24], so the levels of Hurst exponent seen in preterm infants fall within the range of data simulated here. In these simulations, whilst the mean IBI of the external input was equal, the exact distribution of the IBIs was not. However, using the exact same IBI distribution for all simulations, but continuing to vary the Hurst exponent, does not alter the results; the networks driven with LRTCs with higher Hurst exponents exhibit a higher initial rate of change of network parameters (see Fig. S3). This demonstrates that the change in connectivity formation is directly related to the the temporal dynamics of the external input.

Varying the plasticity parameters
In the simulations to this point, LTD of the likelihood of losing/gaining a connection was set to be slightly stronger than long-term potentiation (LTP) (i.e., A D = 0.55 and A P = 0.5, see Methods). This led to the network on average losing connections. We next explored changes in connectivity when these parameters were varied so that potentiation and depression were  Small-world index (average Hurst exponent estimated using DFA in preterm infants is 0.68 [24]).

Varying network size and density
Finally, we investigated changes in the network size and density. Previous simulations have all had 200 neurones starting from a random connectivity structure with a density of 40%.
Varying the network size did not have a major e↵ect on the network parameters, but reduced the variance see Fig. S6. In contrast, varying the initial connection density had a notable e↵ect on the final network structure, Fig 5. With a low initial connection density of 10%, although there is a di↵erence between the network driven with LRTC and that driven by the external input with random burst ordering in terms of the clustering coe cient, the di↵erences in normalised mean path length between the two lead to both networks having similar smallworld indices throughout the simulations. With a high initial connection density of 50% there is a much slower rate of change of the clustering coe cient and small-world index when the network is driven with LRTCs compared with networks with an initial connection density of 40%. Thus, we importantly note that the di↵erence we observe between networks driven with bursts which exhibit LRTCs and those with random temporal ordering is not always observed but is dependent on the network parameters, including the initial connectivity and the plasticity parameters.

Discussion
Here we examined an activity-dependent neuronal network model of connectivity formation in the developing brain. We demonstrate the emergence of small-world topology and the presence of hub neurones, characteristics of brain networks [1,2,29,30]. Furthermore, we found that the temporal ordering of the activity driving the network may be importantwhen the bursty input driving the network had IBIs which exhibited LRTCs this resulted in faster network evolution compared with when the network was driven by the same input randomly reordered in time (within a range of network parameters). Moreover, network evolution dynamics were related to the magnitude of the Hurst exponent, with a faster speed of emergence of small-world properties in networks driven by external input with higher Hurst exponents.

Emergence of small-world topology
Starting from networks with initially random connectivity and allowing network connectivity to evolve using a simple Hebbian STDP rule, we observed the emergence of small-world topology and hub neurones. Small-world connectivity has been shown to arise in simple models due to constraints such as e cient neuronal communication and metabolic costs related to neuronal wiring [31], and hub nodes can arise through mechanisms such as preferential attachment [32]. A number of previous authors have also examined more realistic computational models of brain development, including examining axonal growth in molecular gradients [33] and neurite branching [34]. Van Ooyen and Van Pelt showed that a simple model of connectivity formation, with the growth and retraction of circular dendritic and axonal fields based on neuronal activity, results in an equilibrium in connectivity level after an initial overshoot in connectivity consistent with experimental observations [35]. Meisel and Gross further demonstrated that an activity-dependent model of connectivity formation 'self-organised' to a balanced connectivity level irrespective of the initial level of connectivity [36]. More recently, Damicelli et al. found that a local Hebbian plasticity rule allowed a network to reorganise to a modular structure [37]. Thus, all of these models demonstrate the importance of neuronal activity in connectivity formation. However, whilst it is known that neuronal activity is essential in the developing brain, both before and after birth [12,16], we importantly show that the temporal ordering of activity may also play a role in development.

The role of LRTCs
We find that when driven by burst activity which exhibits LRTCs the network evolves quickly to a state with small-world topology. We speculate that this may be important in development, where quickly transitioning to this type of network will allow for e cient integration and segregation of information in the developing brain. This work was motivated by the observation that inter-burst intervals between bursts of activity in the EEG of preterm infants exhibits LRTCs, with an average estimated Hurst exponent of 0.68 [24]. Evidence suggests that the subplate provides essential input to the cortex during this stage of development [12,13] so we made the assumption that the burst dynamics of the external input were from a region such as the subplate driving our cortical network. However, the LRTCs in burst dynamics are also passed to the network itself, suggesting that other models which exhibit LRTCs in network activity may also quickly evolve their connectivity. For example, this type of burst dynamics with LRTCs in the inter-burst intervals can emerge in a system driven with continuous external input in the region of a critical state [38]. Although our motivation for this work comes from observations in preterm EEG, LRTCs in adults have been observed at all levels of the nervous system including in spike timing [26]. Moreover, the LRTCs in the burst dynamics of our model are in the external input which we consider to be a population of neurones rather than from single neuronal firing.

Future directions -understanding pathology in the developing brain
Premature-born infants display altered brain connectivity, indicated by reduced white matter integrity [39,40] and altered resting state connectivity [41,42], at term-corrected age [39,41] and into childhood [40,42]. Tactile, auditory, visual and noxious stimuli all evoke bursts of activity, observed using EEG, in very preterm infants [8][9][10]. It is plausible that unexpected sensory exposure, which could disrupt the temporal patterning of the ongoing brain dynamics, in the premature period relates to the long-term neurological problems observed in children who have been born very prematurely [43]. Indeed, the number of painful procedures an infant receives in the premature period is correlated with altered brain development, including lower white matter integrity, and lower cognitive ability at school age [40]. Whether sensory stimuli alter LRTCs in the ongoing EEG of preterm infants is an open question, and a possible extension of this work would be to determine whether sensory input disrupts connectivity formation, which may lead to a better understanding of the long-term e↵ects of premature birth.
Long-term depression of hippocampal synapses in mature cultures has been shown to result in weaker synapses followed by selective elimination of very depressed synapses [44].
Moreover, synaptic pruning shares similar molecular pathways with long-term depression (LTD) [22]. This suggests that connectivity formation in the developing brain may be altered through STDP-like mechanisms, which forms the basis of the model we have used here.
In some pathological states, including neurodevelopmental disorders such as autism [45] and schizophrenia [46], brain connectivity is altered, with network architecture which has lower clustering and fewer hub nodes [31]. Both hyper-and hypoconnectivity have been observed in children with autism [47,48], and LTD dysregulation has been identified in mouse models of autism, leading to the suggestion that alterations in synaptic plasticity and pruning may prevent proper development of brain connectivity [22]. Further exploration of STDP models of connectivity formation, such as the one presented here, may shed light on these disorders.

Limitations
A significant limitation of this study is the small e↵ect size that was observed. Moreover, the e↵ect was not observed across all parameter choices, for example, with low initial connection density and with low values of the decay constant of the likelihood of gaining or losing connections. Nevertheless, this work demonstrates the possibility that the temporal patterning of external input can lead to di↵erences in structural connectivity formation in neuronal networks and therefore provides motivation for future empirical and experimental investigations to further elucidate the role that the temporal ordering of activity may play in development.
To simplify the model, all neurones within the network were driven by the same external input. Whilst bursts of activity in the preterm EEG known as delta brushes can occur in a di↵use pattern over large cortical areas, this is not always the case, with localised delta brushes also observed [11]. It would not be realistic to model the whole cortex as being driven by one external input, and a useful extension would be to determine how connectivity changes if neurones were to receive di↵erent external inputs. As well as the external input, neurones that were connected to each other received input when these neighbouring neurones fired.
To avoid either saturation of network firing or quiescence, the weights between individual neurones were evolved according to the level of network connectivity, which can be thought of as a form of homeostatic plasticity [27]. However, a limitation of our approach was that connections could be continually lost (and gained), which can lead to the network forming disconnected components. A possible extension of the work would be to see if maintaining the level of connections within the network, as was done by Damicelli et al. [37] (another form of homeostatic plasticity), but allowing connectivity to continue to evolve under the dynamics of the network, would lead to small-world properties.

Summary
In conclusion, early spontaneous and sensory driven activity is known to be crucial for the development of connections within neural networks. Here we investigated whether the temporal ordering of burst activity within an external driving input a↵ects connectivity formation in a neuronal network model. Using a STDP model of connectivity formation, we observed that the presence of LRTCs in the ordering of burst activity facilitates the emergence of small-world topology and hub neurones. We suggest that early brain activity, driven by the subplate, leads, through activity-dependent mechanisms, to a small-world cortical network structure, and that LRTCs may play an important role in this connectivity formation warranting further investigation.

Neuronal dynamics
Individual neuronal dynamics were modelled as leaky integrate-and-fire neurones described by the di↵erential equation where V is the membrane potential of the neurone, V r is the resting potential, g L is the leak conductance and I(t) is the input (both external input and input from other neurones within the system). When the neurone reaches a threshold membrane potential V thres it fires and is reset to V reset . For all neurones in the simulations V thres = 54 mV , V r = 70 mV , V reset = 60 mV [50]. The leak conductances were randomly chosen from a normal distribution with mean 0.025 and standard deviation 0.005. This heterogeneity in the conductances leads to heterogeneity in the firing dynamics.

Neuronal input
The external input to the system was constructed using fractional Gaussian noise; an example of a process that exhibits LRTCs, with a Gaussian data distribution. The IBI sequence was constructed from a random normal distribution with mean ⇡ 4.5 and standard deviation ⇡ 3.
This sequence was ordered according the ordering of the fractional Gaussian noise process, generating LRTCs in the IBI sequence. The external input to the system was constructed from this sequence of IBIs, with bursts with a duration of 5 and an amplitude of 0.8mV between each IBI, see Fig. 1A. We confirmed that the sequence of IBIs constructed in this way exhibited LRTCs indicated by a DFA exponent greater than 0.5 (see Estimation of the Hurst exponent). The Hurst exponent of the IBIs was altered by varying the exponent of the fractional Gaussian noise.
We compared the connectivity changes, to connectivity changes within a network evolving under external input with random burst occurrence. This input was constructed by randomly shu✏ing the IBIs from the original external input. In this way the two inputs are identical in terms of the distribution of the IBIs (and bursts are identical throughout) and it is only the temporal structure of the input that is altered. We also compared the connectivity changes in networks driven with periodic external input. For this case the IBIs were all identical and were set to a duration of 4. Bursts were identical to the external input which exhibited LRTCs, and had a duration of 5 and an amplitude of 0.8mV In Supplementary Figure 1, the rate of input is calculated within moving windows of length 200 by summing the external input within the window and dividing by the window length.
As all bursts have equal amplitude, this reflects the temporal ordering of the burst activity within the external input.
The external input was the same to all neurones. Each neurone also received, at each timestep, an input from presynaptic neurones that had fired at the previous time-step. Synaptic weights were equal for all connections and were set to 2 pN where N is the number of neurones in the network and p is the proportion of connections in the network. This is calculated at each time step, according to the maximum number of connections which is equal to N (N 1) (networks are directed). This update to the synaptic weights can be thought of as a form of homeostatic plasticity -without it the network either stops firing when connectivity falls (the external input alone is not su cient to make the network fire frequently), or the network starts to fire continually as connectivity is increased. Using this approach, the average levels of activity in the network were maintained across the course of the simulation.

Estimation of the Hurst exponent
The presence of LRTCs in data can be determined through estimation of the Hurst exponent, H. A Hurst exponent of H = 0.5 indicates that the data does not exhibit correlations or exhibits short-range correlations only (for example, white noise). A Hurst exponent of 0.5 < H < 1 indicates that the data exhibits LRTCs. Here we estimated the Hurst exponent using detrended fluctuation analysis (DFA) [51] which calculates the exponent as the slope of the line of best-fit of the average root-mean-square fluctuations across di↵erent box sizes (see Fig. 1C for an example plot and Peng et al. [51] for detailed methodology). Briefly, the signal is first integrated and then divided into boxes of equal length, n. For each box a least-squares fit to the data is found and the integrated signal is detrended by subtracting this local trend.
The root mean square fluctuation, F (n) is calculated and the process is repeated for di↵erent box sizes and the average fluctuation is compared to box size on a double logairthmic plot.
The minimum box size was set to 5 IBIs and the maximum to one tenth of the length of the IBI sequence (the recommended maximum window size [52]). This approach has been used by a number of previous authors to determine the presence of LRTCs in data, including neurophysiological data sets [24,[53][54][55][56][57][58]. As the external input was constructed to have true LRTCs, and these will not be contaminated by noise, it is reasonable to use DFA to calculate a single linear fit to the data. However, in the absence of such prior knowledge, it is more robust to to use maximum likelihood techniques along with model selection methods [59,60].
To examine network firing and periods of activity/inactivity within the network dynamics itself we separated activity using the method of Benayoun et al. [61]. Briefly, two consecutive spikes within a network are separated as distinct bursts if the time di↵erence between them is greater than the average time, dt, between consecutive spikes within the total simulation.
Thus, a single burst consists of consecutive spikes which are less than dt apart. Benayoun et al. used this approach to define avalanches -cascades of network activity -which have periods of separation between them. In this way, avalanches are the same as bursts of activity and so the same approach can be used to determine the bursts here. However, it is worth noting that the term neuronal avalanche is used to define specifically bursts of activity within a network where the distribution of avalanche sizes follows a power-law [62][63][64].
The average IBI within the network was lower than the average IBI within the external input. This meant that, for the same simulation length, there are more IBIs within the network dynamics, so when calculating the DFA exponents (Fig. 1) the maximum box size is larger in the case of the network firing dynamics than for the external input.

Connectivity formation
All networks were initially randomly connected with 40 % connectivity, apart from in section Varying network size and density where this initial connection density was varied. All networks were directed, and connections were also formed and lost in a direction dependent manner. Connections were updated depending on activity within the network, comparing the firing times between all neuronal pairs. Let L be a matrix of values where L(i, j) indicates the likelihood of losing/gaining a connection from neurone i (presynaptic) to neurone j (postsynaptic). L(i, j) was modified by following a spike in neurone j, where t i is the time since the last spike in neurone i, ⌧ is a decay constant and A P > 0 is the amplitude change when t i = 0, and by ◆ following a spike in neurone i, where t j is the time since the last spike in neurone j and A D > 0 is the amplitude change when t j = 0.
A connection from i to j was gained (immediately, if there was not already a connection present) when L(i, j) increased beyond the threshold value g = 2. A connection from i to j was lost when L(i, j) decreased beyond the threshold value l = 2. In order to better take into account temporal dynamics within the system (for example if two neurones only spike together rarely) the values of L decayed with rate ⌧ L . Thus, at each time-step: The loss-likelihood was initially set to zero for all connections and as in Song et al. [28] depression in initial simulations was set to be slightly stronger than potentiation with A P = 0.5, A D = 0.55. We also set ⌧ = 10 and ⌧ L = 100, with relatively slow decay of the likelihood values L allowing for the temporal dynamics of a number of spikes to be taken into account, and the decay of the spike timing ⌧ relatively fast so as to allow the temporal dynamics of the external input to take e↵ect (the periods within the external input are relatively short and so to have an e↵ect the 'memory' within the system must be of a similar level). We consider changes in the results when these parameters are varied in the final section of the Results. Given any two neurones (or, more generally, nodes within a graph) the shortest path length is the shortest distance needed to be traversed to pass from one node to the other. As we set all the synaptic weights as equal, the shortest path length is equivalent to the lowest number of connections between two neurones. The mean path length is then calculated as the average path length for all pairs of neurones within the network [65]. Given two neurones both connected to a third neurone, the clustering coe cient indicates the likelihood that these two neurones are themselves connected [65]. A random network has a low mean path length and a low clustering coe cient [66]. A number of studies have shown that the neural networks have a similar mean path length to a random network (of the same size and density) but are much more clustered indicating that the brain is a small-world network [2,29,31,66,67].
Mean path length and clustering coe cients of the networks were calculated using the Brain Connectivity Toolbox, using the functions for binary directed networks [65].
The values of the mean path length, L, and clustering coe cient, C, are only really meaningful when compared to the average values of random [68,69] or regular [66] networks of the same size (number of connections and number of neurones). We therefore calculated the values for a random network (C rand and L rand ) by averaging over the values for 50 random directed networks of the same size. Random networks were constructed in this way for comparison every 100,000 simulation steps (at the same points as the network clustering coe cient and mean path length were calculated). The clustering coe cient and mean path length of a network from the simulations were then normalised by dividing by the average values from the random networks to obtain the normalised clustering coe cient and normalised mean path length respectively. From these values the small-world index can be calculated [68]: = C/C rand L/L rand (0.2) A value of = 1 indicates that the network is random, whereas > 1 indicates that the network has small-world properties.
We also compared network degree between the simulations by calculating the degree of the neurones, defined as the number of connections that a neurone makes. As the networks were directed we compared the in-degree distribution (the number of presynaptic neurones a neurone has) and the out-degree distribution (the number of postsynaptic neurones a neurone has) separately.