We used a model of interconnected neural masses, where each neural mass represents a large population of connected excitatory and inhibitory neurons generating an EEG or MEG like signal. The model was described in Ponten et al. (2010). The basic unit of the model is a neural mass model (NMM) of the alpha rhythm (Lopes da Silva et al., 1974; Zetterberg et al., 1978). The same model was used in a previous study on bifurcation phenomena of the alpha rhythm (Stam et al., 1999). As previously described, this model considers the average activity in relatively large groups of interacting excitatory and inhibitory neurons. Spatial effects are ignored in this model; we will introduce spatial effects later by coupling several NMMs together. The excitatory and inhibitory populations of each NMM are characterized by their average membrane potentials V_e(t) and V_i(t), and by their pulse densities, i.e., the proportion of cells firing per unit time E(t) and I(t). Static non-linear functions S_E(x) and S_I(x) relate the potentials V_e(t) and V_i(t) to the corresponding pulse densities E(t) and I(t). The excitatory post-synaptic potential (EPSP) and inhibitory post-synaptic potential (IPSP) are modeled by the impulse responses h_e(t) and h_i(t). The constants C₁ and C₂ describe the coupling from excitatory to inhibitory and from inhibitory to excitatory populations respectively. P(t) is the pulse density of an input signal to the excitatory population. Following Zetterberg et al. (1978) the following impulse responses were used:

For h_e(t) the parameter values were: A = 1.6 mV, a = 55 s⁻¹, b = 605 s⁻¹. For h_i(t) the parameter values were: A = 32 mV, a = 27.5 s⁻¹, b = 55 s⁻¹. The sigmoid function relating the average membrane potential, V_m, to the impulse density was also taken from Zetterberg et al. (1978):

Here the parameter values used were: q = 0.34 mV⁻¹, V_d = 7 mV, g = 25 s⁻¹. For the coupling constants we used C₁ = 32 and C₂ = 3 (Lopes da Silva et al., 1974). A schematic representation is shown in Figure 1A. All model parameters are summarized in Table 1. The impulse response and sigmoid functions are shown in Figure 1C.

The final model consisted of several of the NMMs as described above, which were coupled together. Coupling between two NMMs, if present, was always reciprocal, and excitatory. The output E(t) of the main excitatory neurons of one NMM was used as the input for the impulse response h_e(t) of the excitatory neurons of the second NMM; the output E(t) of the second module was coupled to the impulse response h_e(t) of the excitatory neurons of the first NMM. Following Ursino et al. (2007) we used a time delay (T × sample time, with n an integer, 0 < T < 21) and a gain factor. In the present study, n and gain were set to 1 for all connections. A schematic illustration of the coupling between two NMMs is shown in Figure 1B.

In the present study, the model was programmed in Java and implemented in the program BrainWave (version 0.8.47) written by C. J. Stam (Stam et al., 1999). The Java code was based on the Pascal source code described by Schuuring (1988). For the present study the model was extended in order to be able to deal with activity dependent evolution of connection strength between multiple coupled NMMs. The impulse responses, h(t), were implemented as a convolution in the discrete time domain in a similar way as in the Pascal program.

The average membrane potential of the excitatory neurons V_e(t) of each of the NMMs separately was the multichannel output. The sample frequency was 500 Hz. In the present study each run consisted of 9096 (18.19 s). The adjacency matrix at the end of each run was subjected to topographical analysis. Table 1 gives an overview of model parameters and initial settings. These parameters go back to large number of studies with this lumped model, and ultimately to the original model of Lopes da Silva et al. (1974). For the new parameters in the present study, to be discussed below, we choose to keep the parameters of GDP constant (such that they would result in a reasonable outgrowth of connections), and systematically vary the SDP stepsize, as shown in Figures 5 and 6.

In our previous study the strength of connections between neural masses was always symmetrical, and the same strength was used for all connections (Ponten et al., 2010). In the present study we assign a weight w_ij ∈ [0,1] to the connections between two neural masses i and j, where w_ij = w_ji. We update the connection weights w_ij once every 100 time steps using two processes.

The first process for updating the connection weights is GDP. This process reflects that neural masses will increase or decrease their connection strengths to other neural masses when they deviate from a reference value, determined by a distance dependent decay. GDP is also applied once every 100 time steps, directly after SDP. In contrast to SDP, GDP is applied to all connections, even those with w_ij = 0. The connection weight update is given by:

Here, a_GDP is the GDP step size, which is chosen as 0.001. Θ is a modified Heaviside function, with Θ(x) = 1 if x < 0, and Θ(x) = −1 if x > 0. The term w_th = e^−cDist is the reference value to which the weight w_ij is compared, with c determining the decay rate of the exponential. We chose c = 0.2. Dist is the distance between the neural masses i and j, an integer taken from the interval [1,N/2], with N the number of neural masses and circular boundary conditions. The noise parameter η is a real number taken randomly from the interval [0,1]. The GDP and SDP functions are shown in Figure 1C.

The second is SDP. With SDP the weights of all connections w_ij > 0 are updated with a small step given by a Hill function:

Here a_SDP is referred to as the SDP step size and r is the correlation between the pulse densities E(t) of two neural masses i and j, computed over the preceding 20 time steps, +1. The range of r is [0,2]. The exponent b, which determines the steepness of the Hill function, is chosen as 2. H determines where the function crosses the line Δw_ij = 0 and is chosen as 1. With this function, Δw_ij is bounded between −0.5 a_SDP and 0.3 a_SDP. SDP reflects how the connection strength between two neural masses will increase if they show correlated firing, and will decrease if they do not. If the update results in w_ij < 0, w_ij is set to 0, and the connection is considered to be lost. If the update results in w_ij > 1, it is reset to w_ij = 1.

The weighted adjacency matrices of the model were analyzed with weighted graph analysis. Each neural mass was considered as a node, and the weight matrix with w_ij = w_ji determined an undirected weighted graph. For each node, the weighted clustering coefficient C_w,i was calculated following Stam et al. (2009) as:

The average weighted clustering coefficient was determined by averaging C_w,i over all N nodes. The weighted shortest path was determined by the harmonic mean of the shortest paths.

The length of a path l_ij between nodes i,j was defined as the sum of the inverse of the weights of all edges making up the path. In the case of w_ij = 0, we took 1/l_ij = 0. The weighted clustering coefficient and path length were both compared to the average clustering coefficient and path length of 50 surrogate networks, preserving the number of nodes and edges, the degree and symmetry of the graph, to obtain the normalized measures gamma and lambda.

The weighted degree correlation or assortativity coefficient was based upon the formula of Leung and Chau (2007):

where H is the total weight of all links in a network. In this formula, ϖ_φ is the weight of the ϕth link, F(ϕ) is the set of two vertices connected by the ϕth link, and k_i is the degree of vertex i. The weighted assortativity coefficient R_w scales between −1 and 1.

The modularity of the weighted connection matrix was determined using a modification of the approach by Guimera and Nunes Amaral (2005), adapted for weighted networks. The weighted modularity index Qmw is defined as:

Here, m is the number of modules, l_s is sum of the weights of all links in module s, L is the total sum of all weights in the network, d_s is the sum of the strength of all vertices in module s. A simulated annealing algorithm was used to find the optimal way to divide the network into modules. Initially, each of the N nodes was randomly assigned to one of m possible clusters, where m was taken as the square of N. At each step, one of the nodes was chosen at random, and assigned a different random module number from the interval [1,N]. Modularity Qmw was calculated before and after this. The cost C was defined as −Qmw. The new partitioning was preserved with probability p:

Here, C_f is final cost and C_i is initial cost. The temperature T was 1 initially, and was lowered once every 100 steps as follows: T_new = 0.995 T_old. In total, the simulated annealing algorithm was run for 10⁶ steps.

In the first experiment the influence of the GDP in the absence of any SDP was studied for networks with size N = 32 and 64. The results are shown in Figure 2. In Figure 2A the evolution of gamma, lambda, degree correlation, and weighted modularity of the adjacency matrix are shown for the first 100 epochs (each epoch consisted of 9096 time samples), averaged over 10 runs (error bars ± 2 SD). After an initial transient, the normalized clustering coefficient gamma rises to a stable value slightly above 1.2. The normalized path length lambda drops slightly in the beginning and then stays close to 1. The degree correlation evolves from values close to 0 to values that fluctuate around 0.1 (please note that in Figure 2 the scale for R_w on the y-axis is 2 × [R_w + 1]). Weighted modularity starts at 0, and evolves toward a stable value around 0.236 (value for corresponding random network: 0.091). The final weighted adjacency matrix is shown in Figure 2B. The exponential decrease of connection strength as a function of distance can be clearly seen. Although the connection strengths of the different neural masses are slightly different due to the random adjustment steps of the connection strength, this is not visible in Figure 2B. The networks in Figures 2B,D are not really lattices, even though they look quite regular. Please note that the error term in Eq. 3 introduces a small but important amount of noise in the update of the connection weights. The network in Figure 2B is in fact closer to a small-world network, in agreement with the high gamma and low lambda.

The final node strength (sum of all connection strengths) distribution P(S) is shown in the inset in the right upper corner of Figure 2. This confirms that all nodes have about the same strength. Figures 2C,D show the results for the larger network of N = 64. Note that the formula for the GDP depends upon absolute, not relative distance; we used this approach to keep the average node strength within certain limits, irrespective of network size. Results for the large network are similar to those for the small network with a few exceptions. The normalized clustering coefficient rises to substantially higher values than in the small network, and reaches values around 2.2. The degree correlation fluctuates around 0, and does not take on positive values as was the case for the small network. The weighted modularity reaches a value of 0.403 (random network: 0.149).

In conclusion, starting from an initially unconnected network, GDP gives rise to a network with connection strengths that drop exponentially with distance, have typical small-world features but a narrow strength distribution and no visual indication of modular structure, although values of weighted modularity are higher than those of random networks.

In this experiment the effect of SDP was studied in isolation. Since SDP can only change the strength of existing connections, and does not involve the growth of new connections, an unconnected network could not be used as the initial state. Instead we used a network with average degree k = N/2, and a fully random distribution of connections over the network. All connections weights were random numbers from the interval [0,1]. For the SDP step size parameter, a_SDP, a value of 0.005 was chosen. In Figure 3A the plots for N = 32 are shown. While lambda stays around 1, gamma evolves toward a stable value of 1.2. The degree correlation rises to a positive value around 0.1. Weighted modularity reaches a value of 0.229 (random network: 0.175). The strength distribution (inset right upper corner) is rather broad. The final weighted connection matrix, shown in Figure 3B, reveals a pattern of widely scattered small clusters of connection strengths close to 1, surrounded by a background where the connection strength is close to 0 (please note that there has been no reordering of the nodes in this or any other of the figures). Results for the large network are shown in Figures 3C,D. While the general pattern is the same as for the small network, there are a few differences. Gamma and the degree correlation evolve to higher values, and still are increasing after 50 epochs. Weighted modularity reaches a value of 0.354 (random network: 0.161).

Thus, SDP operating on a network with initial random connection topology and strength, results in a small-world network with scattered clustering, assortativity and broad strength distributions. The latter features were either absent or inconsistent with networks based upon GDP alone. SDP alone did however not result in networks with clear modules, although the value of weighted modularity was higher than for corresponding random networks.

The central experiment involved the combined effect of GDP and SDP, starting from an initially unconnected network. In these experiments GDP step size was fixed at 0.001, while SDP step size varied from 0 to 0.012. First, in Figure 4 we show the evolution of small and large networks for an SDP step size of 0.008, and 100 epochs (100 × 9096 time steps), averaged over 10 realizations. Figure 4A shows that lambda stays around 1, except for a brief small increase around epoch 20. Gamma is initially high, drops very quickly, and then increases again around epoch 20, first quickly, then slower, to reach a value of 1.4. The degree correlation evolves toward positive values, and has a temporary maximum around epoch 20. Weighted modularity increases up to 0.278 (random network: 0.109). The final strength distribution is slightly broad. The final weighted connection matrix, shown in Figure 4B is of special interest. We can clearly see a complex pattern with a dense module in the upper left corner, a smaller overlapping module in the middle, and a large but fragmented module occupying the lower right corner. Results for N = 64 are shown in Figures 4C,D. The pattern is the same, but higher final values are obtained for gamma (2.37), weighted modularity (0.450; random network: 0.147) and the degree correlation (0.113). Also, the peak in the degree correlation around epoch 20 is more outspoken than in the small network. The weighted connection matrix (Figure 4D) shows a delicate pattern of multiple, partially overlapping modules. The appearance of the weighted connection matrix for small and large networks for different values of the SDP step size is illustrated in more detail in Figure 5. The transition of a mostly distance dependent connection pattern for SDP step size = 0.002 (upper and lower leftmost panels) toward a more complex modular pattern for SDP step size = 0.010 (upper and lower rightmost panels) is clearly visible.

The weighted modularity index Qmw was also studied for the small network, 200 epochs, and SDP stepsize varying from 0 to 0.012 in steps of 0.001. The average result (error bars ± 2 SD) of 10 networks is shown in Figure 6. Qmw shows two maxima as a function of SDP step size, one for step sizes of 0.003–0.004, and higher one for step sizes of 0.009–0.010. Of interest, the first maximum corresponds to a division of the network into three modules, while the larger second maximum corresponds to a division of the network into two modules. The lower curve shows the results for random networks with the same average connection strengths. Except for SDP stepsize 0.004 and 0.005, and SDP stepsize > 0.010, the error bars (error bars ± 2 SD) do not overlap.

Network development with an appropriate balance between GDP and SDP thus gives rise to small-world networks with high values of kappa, assortative degree coupling, and a complex modular structure. Size and number of modules depend upon the balance between GDP and SDP, as well as network size.

An important question is how plastic brain networks recover from a sudden brain lesion. Lesions were modeled by setting all connections of neural masses 1–5 to very low values (0.1 × η, with η a random number from the interval [0,1]) at epoch number 51. Network parameters and evolution for epoch 1–50 were the same as in Experiment 3. Figure 7A shows the results for the small network. The first part up to epoch 50 shows a similar evolution as Figure 4A. The lesion at epoch 51 gives rise to a sharp increase in lambda, gamma, and the degree correlation; weighted modularity shows a modest decrease. During the next 50 epochs lambda, gamma, and modularity recover to their original values. The degree correlation shows a temporary increase after the lesion followed by a slow recovery of the values before the lesion. The final weighted connection matrix is shown in Figure 7B. It displays a complex pattern, with multiple partially overlapping modules of different size. The results for the large network, shown in Figures 7C,D, agree with those for the small network. A more detailed view of connectivity changes induced by the lesion is shown in Figure 8. The upper left corner shows a network of N = 32 after an evolution of 100 epochs, with the same settings as before. The upper middle pattern shows the immediate lesion effect: all connections of node 1–5 are weakened to very low values. The subsequent panels (upper right, lower left to right) show the network changes in steps of five epochs. There is a clear outgrowth of new connections, resulting in a partial restoration of the module in the upper left corner.

The first process we included in our model was that of GDP. There is increasing evidence that neurons form highly dynamic elements, not only in terms of their electrical activity, but also in terms of their ever changing morphology. It has been suggested that neurons may regulate the outgrowth or retraction of dendrites and axons in such a way that they maintain a certain level of activity (Butz et al., 2009a). Such homeostatic aspects of structural plasticity could also determine the response of networks to lesions (Butz et al., 2009b). In our model we incorporated this in an indirect way, by evolving connections strengths such that the mean strength (sum of all connection weights of a node) of each node was kept within certain limits. The mean connection strength is a key determinant of the activity level of neural models, and limiting this strength therefore ensures that the activity level remained within certain limits.

A second important consideration is space. In many graph models space is not considered explicitly. However, the topology of developing and mature brain networks is hard to understand without taking distance into account. Kaiser and Hilgetag (2004) have shown that random, distance dependent outgrowth of connections in a two-dimensional geometric model can explain many of the topological features of real anatomical brain networks of cats and macaques. However, to explain the emergence of modules, ad hoc assumptions about the timing of connection growth had to be made (Kaiser and Hilgetag, 2007). More recently it has been shown how the exponential distribution of connections in brain networks can be explained by a random outgrowth of dendrites and axons (Kaiser et al., 2009). We incorporated this in the exponential shape of the GDP function (Eq. 3). The randomness in the process was reflected by the small, random update steps of GDP. As a result, the connections strengths in our model showed small, random, fluctuations around an average exponential decay as a function of distance (Figures 2B,D). GDP alone can explain how an initially unconnected network can evolve toward a network with small-world features, in agreement with the observations of Kaiser and Hilgetag (2004). However, GDP does not consistently produce assortative networks, and does not fully explain modularity. The presence of relatively high modularity in non-modular lattice-like networks is a well-documented conceptual shortcoming of the modularity metric. Although the final values of the weighted modularity were higher than those of a random network, no modular structure was visible in the connection matrix. Also, GDP only depends upon the connection topology and weights, and does not take the activity on the network into account.

In contrast to GDP, SDP represents a true interaction between dynamics on the network and the underlying structural connection strengths. At the neuronal level, much has been learned about such mechanisms (Abbott and Nelson, 2000). Depending upon their correlated activity levels the synaptic strength between neurons can be adapted by pre- and post-synaptic mechanisms such as synaptic scaling and STDP. In addition, rewiring of existing connections and outgrowth of dendrites and axons play a role at the neuronal level (Butz et al., 2009a). Thus, both functional as well as structural mechanisms play a role with synaptic plasticity. A central idea is that of Hebbian learning: connection strengths between neurons that show correlated firing increase over time. While Hebbian learning may explain selective reinforcement of functionally important connections, this mechanism may be inherently unstable. Instability could be dealt with by synaptic scaling, and depletion of pre-synaptic resources (Levina et al., 2007). The mechanism of SDTP combines selective reinforcement with stability since excessive connectivity will disrupt the exact timing of pre- and post-synaptic neurons (Song et al., 2000). Models incorporating one or more of the synaptic mechanisms of plasticity typically evolve from an initial random topology toward a complex architecture with typical small-world features (Siri et al., 2007). In addition, SDTP and pre-synaptic depletion may explain how such evolving networks become critical (Shin and Kim, 2006; Levina et al., 2007; Sendina-Nadal et al., 2008). However, it is not clear how these findings should be translated to the macroscopic level of interacting cortical regions and observed correlations in EEG, MEG, and fMRI studies. Also, models at the neuronal level do not seem to explain the emergence of modular structure and assortative mixing.

The SDP in our model is an attempt to describe how the various mechanisms of functional and structural plasticity described above might be reflected at the level of interacting neural masses. We assume that neural masses will increase the strength of their excitatory interconnections when they oscillate more synchronously, and decrease this strength otherwise. This idea was implemented in the Hill type SDP function (Eq. 4). At this stage we cannot derive the SDP function from the multitude of synaptic mechanisms involved. However, it seems reasonable to assume that increased correlated firing at the neural level will be reflected by increased synchronization at the population level. We used the spike densities, not the average voltages, to determine the correlated activity of two neural masses, since the spikes may be closer to the actual neural interaction processes. SDP cannot give rise to new connections, so if we want to study its effect we have to assume some pre-existing topology and distribution of connection weights. This was explored in Experiment 2, where we applied SDP to sparsely connected networks with random topology and random connection strengths. SDP resulted in a decrease of some connections and an increase of other connections. Interestingly, SDP resulted in scattered small clusters (Figures 3B,D). SDP is a positive feedback process, and it seems likely that small groups of highly synchronous and connected nodes will tend to grow and recruit other nearby nodes in the process. This might explain why the resulting networks are small-world, but also assortative (Figures 3A,C). However, SDP does not result in true, large-scale modules, even though the weighted modularity is higher than that for a corresponding random network. In addition, this SDP-induced modularity is probably neurobiologically unrealistic, as the modules in Figure 3 seem to be formed between “spatially remote” regions. Furthermore SDP is an inherently unstable mechanism, which drives all connection weights to either their lowest allowed value (0) or their highest possible value (1). The instability of SDP without GDP is also clear from the fact that gamma and R_w do not reach stable values even after 100 runs (Figures 3A,C). Therefore, SDP alone does not provide a realistic account of mechanisms of plasticity, and it does not account for the modularity found in real networks.

The central idea of our model is that we need both GDP and SDP to simulate how activity and connectivity interact in developing plastic neural networks. This interaction is demonstrated in Experiment 3. Varying the strength of SDP for a constant strength of GDP gives rise to complex networks with small-world features and, most importantly, assortative mixing in combination with clear modularity (Figure 4). Modularity for the network with combined GDP and SDP was higher (0.278) than that for either GDP (0.236) or SDP (0.229) in isolation. Only the combination of GDP and SDP resulted in a clearly visible modular structure. The fact that modular networks emerged for a range of SDP strengths can be understood by considering how the two processes interact. For absent or very weak SDP the process is obviously dominated by GDP alone, resulting in non-modular networks (Figure 5). When SDP is too strong, any new weak connections generated by the GDP mechanism will be removed by the SDP function due to a weak correlation between firing patterns. Only for intermediate values of SDP (between 0.002 and 0.010) we have a situation where some of the new connections generated by GDP are strong enough to cause synchronization of neural masses, and subsequent strengthening by the SDP function. Interestingly, modularity as a function of SDP step size has two maxima, one at 0.003–0.004 and a second one at 0.009–0.010. It seems the dynamics of GDP and SDP interact at multiple “stable states”, each corresponding to a different number of large-scale modules.

Only a few other studies have attempted to explain the emergence of modular structure in complex networks. Fuchs et al. (2009) investigated synchronization dynamics in modular networks, but did not consider how the modular topology could arise in the first place. Aoki and Aoyagi (2009) studied the interaction between network topology and synchronization dynamics in a network of coupled phase oscillators. A two-cluster state could emerge in this model, suggestion that connectivity/dynamics interactions may be relevant for the emergence of modularity. However, the details of this simple model are hard to translate into possible neurobiological mechanisms of module formation. In a series of studies it has been shown that synchronization dependent rewiring of networks of coupled logistic functions can give rise to modular structure (van den Berg and van Leeuwen, 2004; Zhou et al., 2007; Rubinov et al., 2009b). While this is an elegant minimal scenario, it requires chaotic dynamics of the units, which may not be a realistic scenario for the dynamics of neural masses. We have previously shown that the NMMs used in our study have dynamics that is at most weakly non-linear, and not chaotic (Stam et al., 1999). In addition, the rewiring scheme that was used in the previously mentioned studies is biologically unrealistic since it does not take into account the influence of distance on neural connectivity, and does not handle outgrowth of new connections and loss of old connections. The model of Kaiser and Hilgetag (2007) does take into account the importance of distance, but, in order to explain modules, required the extra assumption of specific time intervals during which connections in different parts of the system grew. In our model, modularity evolved naturally without such timing assumptions.

A combination of GDP and SDP gives rise to complex modular networks. Even in the stable state, with a relatively constant modular architecture, these networks are still dynamic, with the slow appearance and disappearance of connections and the constant adjustment of connection weights. It would be interesting to know how this ongoing plasticity would affect the response of the brain network to lesions, to better understand observations on network characteristics in neurological patients (Butz et al., 2009b). We have shown that an acute “lesion”, modeled by a sudden loss of connections of a few adjacent nodes, results in an immediate change in network topology with increased path length and clustering, increased degree correlation and decreased modularity. Over time, the network recovers most of its original structure, including the modularity, but the recovery rate and pattern is different for different network properties. While gamma, lambda, and modularity show an exponential approximation to the original values, the degree correlation shows a transient positive peak some time after the lesions. This peak is reminiscent of the peak that is also observed in the initial development, at the stage where the modular structure is rapidly evolving. Recovery from a lesion thus reflects to some extent a replay of events during network evolution. Other studies that have investigated the effect of lesions on complex networks have mainly concentrated on acute effects, and have not yet taken into account the role of plasticity in network recovery (Honey and Sporns, 2008; Alstott et al., 2009; Rubinov et al., 2009a). The study by Butz et al. (2009a,b) did take into account plasticity after lesions, and also reported a slow recovery of the initial network. However, this study involved networks of neurons, and does not describe changes in network topology on larger scales.

Modeling of lesion effects in plastic networks with activity–connectivity interactions could help to understand empirical observations with EEG, MEG, and fMRI in patients. There is now increasing evidence that brain lesions, in particular brain tumors but also traumatic brain injury, may have widespread effects on networks, and may disrupt the normal network topology (Bartolomei et al., 2006a,b; Bosma et al., 2009; Nakamura et al., 2009). The nature of these changes may depend upon time of observation and the frequency band studied, as well as therapeutic interventions (Douw et al., 2008). Network randomization has mainly been found in higher frequency bands (Bartolomei et al., 2006b; Bosma et al., 2009). In lower frequency bands, brain tumors may give rise to abnormal increase in clustering (Bosma et al., 2009). Of interest, the latter observation corresponds to the increase in the normalized clustering coefficient observed after lesions in our model.

Obviously, the current model is only a crude initial attempt to understand how activity–connectivity interactions can give rise to complex networks with modular and assortative structure, and how such models respond dynamically to damage. However, a major advantage is that the model makes predictions at the level of macroscopic measures such as EEG, MEG, and fMRI BOLD signals. This creates the opportunity to study fundamental mechanisms in the model, and make predictions with respect to output signals that can be directly tested in healthy subjects and patients. In a previous study we have shown that a consistent but complex relation exists between the topological properties of the adjacency matrix and the topological properties of functional networks derived from the output signals of the model (Ponten et al., 2010). Ideally, an iteration of modeling and empirical studies will increase our understanding of fundamental brain mechanisms underlying network development and plasticity. Major challenges for future studies are linking the evolution functions, in particular the SDP, more directly to known synaptic mechanisms, devising ways to test information processing in the model, and implementation of more sophisticated disease models of growing lesions, white matter disturbances and neurodegenerative disorders such as Alzheimer's disease (Stam et al., 2009).