Regarding the choice of our model, we consider a random graph [24] network based on McCulloch-Pitts neurons as favorable for the examination of rTMS effects on the band power, because we aim at modelling the influence of rTMS on the brain area directly affected by the magnetic field of the TMS coil. The random graph model for the generation of the graph is the G(n,p) model proposed by [24], which is not differing from the exponential graph models (ERGM) if the edges metric constitutes the only network metric [17] . We also prefer the simple McCulloch-Pitts neurons, since a complex model for the behaviour of the neurons based on differential equations like the Hodgkin-Huxley model [25] often exceeds, beyond numerical problems, the limits of practical computability.

The network size has no statistical relevance. This scale-freeness of the model has important consequences: a) the results can be transferred on networks of different sizes and b) we can predict local field potentials, cortical potentials and even potentials on the surface of the head. This complies with experimental findings that indicate the independence of scale [26], i.e. microelectrode recordings of the activity of a small number of neurons and scalp electrode recordings reflecting the activity of millions of neurons are essentially identical.

The results of the simulation itself aim at identifying how band power depends on rTMS frequency and intensity, and, in doing so, advance the development of rTMS as an efficient therapeutic treatment for mental disorders [27] and diseases like tinnitus [28].

The results replicate in vivo experimental observations in three main ways: gamma band power is increased, alpha band power depends on rTMS intensity, and band power is independent of rTMS frequency in most protocols. The model also reveals the stability of the band power in the beta band.

We observe an increase in the gamma band power for all stimulus intensities and frequencies above 0.5 Hz. These findings agree with in vivo experimental results [13], [14] and imply that short rTMS pulses at repetition frequencies greater than 0.5 Hz significantly enhance gamma band activity. No experimental results have been performed with a 0.5 Hz stimulus protocol, and thus this cut-off frequency can neither be confirmed nor substantiated by experiments.

In the alpha band, the simulation shows that changes strongly depend on rTMS intensity. This observation accounts for contradictory experimental results, regarding the effect of rTMS on the power in the EEG alpha band [8], [10], [11].

As seen in recent research studies and confirmed by these findings, the rTMS intensity does significantly affect the EEG power across all frequency bands. The stimulation intensity influences the network activity by causing a higher number of spikes at the neuronal level and accordingly changes of the eLFP with increasing intensity. It is important to note that the effects are studied on a neuronal network level which means that our model neglects potential side effects of the applied rTMS energy.

This model does not explain all experimentally observed results of rTMS. Our simulation and most in vivo experiments do not reveal information about long-term effects of rTMS on the neuronal level e.g. changes of synapses like synaptic plasticity or synchronization of brain activity. This point must be investigated further to predict the efficiency of rTMS therapies.

In addition, we assume a uniform perturbation of the network by the rTMS stimulation. It is left open how the stimulation affects remote brain areas not directly affected by the magnetic field. The effects of rTMS on remote, not directly stimulated cortical areas are not included in this model. We believe that first and foremost, it is vital to understand the effects of rTMS on the directly stimulated brain areas and to close the gap between model behaviour and EEG measurement data.

The discrete time evolution for each neuron is given by the input function I_i(t) and the output function E_i(t) with a threshold value of Θ_i=−60 mV for the generation of an action potential (AP), and with constant c denoting the change in postsynaptic potential triggered by the arrival of an AP at the synaptic cleft [22]. Firing neurons either cause a hyperpolarisation or depolarization of ± (0.2 mV–1 mV) [22]. For sake of simplicity, we assume a value of ±0.5 mV for the effect of inhibitory and excitatory postsynaptic potentials on the cell membrane and set c=0.5. The effects of a single AP on the postsynaptic membrane potential is not strong enough to generate an AP since the membrane potential returns to its resting potential of −70 mV after a few milliseconds [22]. The net effect of spatial summation depends on various parameters [22]. In the model, we assume a linear spatial and temporal summation.

Whether the synapses are of the inhibitory or the excitatory type is coded in the randomly generated adjacency matrix a_ij [24]. The number of synapses depends on the connection probability in the network which in turn determines the average node degree k of each neuron in the network [24]. Existing experimental results could be explained best for k=120. Table 2 shows the connection probability p resulting in k=120 for different network sizes. These theoretical considerations are substantiated by microscopic analysis of mammal cortical tissue indicating a connection probability of p=0.12 for a cluster of 1,000 neurons and a ratio of 8020 between excitatory and inhibitory synapses [23].

For determining the stimulus intensities, we assume that a single TMS pulse triggers the simultaneous depolarization of neurons [30]. Numerical analysis of our network revealed a mean change of 1.25 mV of the synaptic potential per neuron for k=120. Thus, we set the TMS induced disturbance to 160%, 240%, 320% and 400% of the mean change of the synaptic potential which corresponds to a change of 2 mV, 3 mV, 4 mV, 5 mV. One can think of that as the increase caused by the simultaneous spiking of either 4 neighbours (TMS-INTENSITY 4), 6 neighbours (TMS-INTENSITY 6), 8 neighbours (TMS-INTENSITY 8) or 10 neighbours (TMS-INTENSITY 10) [30]. These values guarantee stable oscillatory behaviour of the network after the TMS stimulation for all network sizes. The TMS stimulation itself is applied within a period of 7.99 s, beginning at 8,000 s and lasting to 15,999 s. The simulation output is defined as the sum of the output E_i(t) of all neurons at time t which corresponds to macroscopic properties like the LFP or the surface EEG and will be called equivalent local field potential (eLFP). For the processing of the output ρ(t) (Fig. 4) of the model, we selected one segment of 1 s ending immediately before the first (PRE) and one segment of 1 s beginning after the last stimulation (POST). For all combinations of the model parameters TMS-INTENSITY and TMS-FREQUENCY we calculated 20 independent runs, each of the runs incorporating a new, random adjacency matrix (this lead to 4 * 6=24 models times 20 independent runs which leads to 480 model runs). In rare cases ρ(t) showed limit cycle behaviour for certain model runs, which were treated as missing values. In this case the independent run was excluded from the statistics. This left us with 17 independent runs as a source of variance of ρ(t). As the variance originates from the adjacency matrix, which reflects individual interconnection of the neurons, we argue that this would also reflect variations of the network likely to be found in different cortical columns within the same individual or even across individuals.