Circuits in the motor cortex explain oscillatory responses to transcranial magnetic stimulation

The activity of motor cortical neurons directly contributes to the generation and control of movement in humans. Noninvasive brain stimulation techniques, notably transcranial magnetic stimulation (TMS) have been used to study the connection from cortex to muscle in recent decades (Ni, Müller-Dahlhaus, Chen, & Ziemann, 2011b; Volz, Hamada, Rothwell, & Grefkes, 2015). TMS is also used as a measurement tool for cortical excitability with direct applications in epilepsy and stroke assessment, as well as a potential neuromodulation therapy in depression (Badawy, Strigaro, & Cantello, 2014; Smith & Stinear, 2016). Though there are many examples of clinical experimental research (Blicher, Jakobsen, Andersen, & Nielsen, 2009; Carpenter et al., 2012; Kolbinger, Höflich, Hufnagel, Müller, & Kasper, 1995; Stinear et al., 2007), there currently lacks a mechanistic neurobiological understanding of the interaction of TMS and neuronal circuits.

A single-pulse TMS at a sufficient intensity applied over the motor cortex evokes electrical currents in the neurons and a subsequent contraction in downstream muscles. Following the application of TMS over the motor cortex, high-frequency volleys of electrical activity of 600–700 Hz have been recorded from the epidural space in the cervical spinal cord of humans, as well as cats and monkeys (Di Lazzaro et al., 2012; Ziemann & Rothwell, 2000). These cortically evoked responses in the pyramidal tracts were first described in Patton and Amassian (1954)’s study; however, the exact physiological mechanisms involved in the generation of these waves remain unclear.

Cortically evoked volleys of electrical activity have been shown to be due to the direct and indirect activation of corticomotor cells. The initial volley, known as the D (direct)-wave, is considered to be due to the direct activation of the layer 5 pyramidal tract neurons. Following the D-wave, a series of three to four I (indirect)-waves have been attributed to the activation of upstream cortical neurons with differing synaptic delays in recurrent or polysynaptic connections with the pyramidal neurons (Di Lazzaro et al., 2012; Rusu, Murakami, Ziemann, & Triesch, 2014). I-waves originate from motor cortex circuits and are susceptible to cortical injuries, unlike the D-wave, which is relatively independent of the functional state of the cortex (Patton & Amassian, 1954; Zewdie & Kirton, 2016). Hypothetical models have been proposed to explain I-wave generation involving neuronal circuits of excitatory and inhibitory interneurons, but these have not been explicitly tested in population-based models (Di Lazzaro & Ziemann, 2013; Ziemann & Rothwell, 2000).

The motor cortex is composed of a large number of neurons with different size, location, orientation, and functions (Young, Collins, & Kaas, 2013). The neurons of the motor cortex are structured in a highly stereotyped manner of layers and columns (Douglas & Martin, 2004; Hatsopoulos, Xu, & Amit, 2007). Canonical cortical circuits of excitatory and inhibitory neurons have previously been suggested by Gilbert and Wiesel (1983) and Douglas and Martin (2004), involving recurrent connections between the six layers of the cortex. The circuits in the motor cortex and how they are activated by TMS is difficult to establish through anatomical and electrophysiological measurement alone because it is currently infeasible to record from and manipulate precisely the large numbers of neurons involved. Though there have been some recent efforts to experimentally explore the cellular mechanisms of TMS, computational modelling can test the plausibility of proposed mechanisms and shed light on hypotheses involving complex interactions between larger scale neuron groups (Epstein, 2008; Murphy, Palmer, Nyffeler, Müri, & Larkum, 2016).

The generation of I-waves has previously been modelled using compartmental models of individual neurons and spiking neural networks. Rusu et al. (2014) modelled a population of layer 2 and 3 neurons projecting onto a detailed layer 5 pyramidal cell model and suggested that I-waves are generated by the delay due to synaptic inputs arriving at different parts of the layer 5 dendritic tree. This was a feedforward model and ignored any recurrent connectivity that exists within and between cortical layers (Weiler, Wood, Yu, Solla, & Shepherd, 2008). Esser, Hill, and Tononi (2005) used a large-scale spiking neural network with detailed anatomical and physiological constraints to model the motor cortex and the effects of GABA enhancement, as well as single and paired-pulse TMS. However, this study has not been replicated and did not identify or examine the key contributions of the components of the circuit, rather it treated the whole neural network as a black box without showing which connections or features were important in capturing the TMS response.

Activation of TMS-induced descending volleys is dependent on the orientation of the TMS coil due to the direction of induced current (Ni et al., 2011a). A current in the posterior-anterior (PA) direction preferentially activates early I-waves, which increase in amplitude and number as stimulation intensity increases. Stimulation in the anterior-posterior (AP) current direction results in fewer, lower amplitude waves, which have a longer onset latency. Multi-scale computational modelling exploring the effects of TMS on individual neurons suggests that TMS activates neuronal axons aligned to the direction of the induced electric field, depending on the cortical morphology (Aberra, Wang, Grill, & Peterchev, 2020; Seo, Schaworonkow, Jun, & Triesch, 2016). Shifting the direction of TMS from an AP to PA direction potentially activates different circuits of the motor cortex. The later onset of AP–stimulation induced I-waves is thought to be the result of polysynaptic pathways through the premotor cortex (Spampinato, 2020).

While previous computational models have captured several properties of I-waves, a coherent conceptual framework for understanding the response of cortical circuits to TMS is lacking. The aim of this study was to develop a structural model of the primary motor cortical circuitry and determine whether recurrent connections in the cortex explain physiological responses to TMS as observed by the generation of I-waves in the spinal column. We used a spiking neural network model of excitatory and inhibitory neurons in the different cortical layers, connected with plausible connection probabilities and synaptic strengths. This motor cortex model is described in Haggie, Besier, and McMorland (2023a) and uses elements from a previously published cortical model by Potjans and Diesmann (2014) and the model put forth by Esser et al. (2005).

The model was altered in several ways to explore the contributions of the structural arrangement and functional connectivity between neuron groups to the response of TMS over a region of motor cortex, including testing several circuits previously hypothesised by Di Lazzaro et al. (2012) to be required for I-wave generation. Our model replicates the effects of single-pulse PA and AP TMS, increased GABAergic inhibition, and facilitatory paired-pulse TMS and also shows adaptation to TMS strength and a cortical silent period. To the best of our knowledge, this is the first study to use a large-scale population model of the response to TMS in the AP and PA current directions explaining circuit-level contributions.

A consistent connectivity was defined in the TMS simulation studies by reading in tables of pre- and postsynaptic neuron indices corresponding to a single instance of the cortical model with local connectivity. A single instance of the model was used in most simulations to eliminate the possibility that differences observed were due to specific definition of the randomly chosen connections of individual neurons within the network structure. We observed similar results in the network firing rates and patterns as well as observed responses to a simulated TMS perturbation, even with slight changes to the network structure based on the randomness of input, spatial and connectivity definition in this model.

External input to the model representing activity received by motor cortex neurons from other brain regions was a Poisson-distributed spike train of 8 Hz. The number of these inputs for each neuron was constant at 2,000 inputs for excitatory inputs and 1,850 for inhibitory inputs, as previously used in the layer-independent protocol from (Potjans & Diesmann, 2014). Though thalamic input mostly targets superficial layers 2/3 to 5, inputs from other cortical areas also contribute (Hooks et al., 2013; Okoro et al., 2022). Connections such as from the orbital cortex are shown to target layer 6 (Hooks et al., 2013), therefore a homogenous input was used in this model.

The spontaneous activity of the model, in response to poisson input, shows a resting-state population firing rate of 3.2 Hz in layer 2/3E and 10.5 Hz in layer 5E neuron groups of the cortex, which is similar to intracortical recordings (Dabrowska et al., 2021). Figure 3 shows the spontaneous activity of the layer 2/3E and layer 5E neuron populations in which bands of oscillatory waves of activity can be observed similar to those reported in Esser et al. (2005).

TMS has been shown to induce the activation of single neurons in areas of the cortex measuring less than 2 mm in diameter (Romero, Davare, Armendariz, & Janssen, 2019). The effect of TMS is also thought to act superficially on the cortex with reduced strength at further depths (Rudiak & Marg, 1994; Trillenberg et al., 2012). A single TMS pulse was modelled as a direct input current to a proportion of individual neurons, that is, a value of 25% of fibre terminals stimulated, which is the proportion of fibres stimulated in nonhuman primate studies to elicit motor evoked potentials (Romero et al., 2019). The strength of the TMS activation at different depths was modelled with a linear reduction in the input current to neurons as the depth (Z) increases (Figure 4). Different amplitudes of TMS were modelled by altering the proportion of neurons activated.

A group of neurons representing the premotor groups was added to the cortical model to explore the different pathways activated during AP and PA TMS protocols. The majority of neurons from the premotor area projects to layer 2/3 in the motor cortex with a small proportion of connections to layer 5 (Ghosh, Brinkman, & Porter, 1987; Leichnetz, 1986; Lu, Preston, & Strick, 1994; Matsumura & Kubota, 1979). The response to PA TMS was modelled by the activation of the primary motor cortical neurons, while the AP response to TMS was modelled by stimulation of the premotor neuron group, shown in Figure 5. Previous computational models have suggested that TMS activates axonal collaterals aligned to local electric field direction and that AP stimulation results in activation of collaterals directed toward the premotor cortex, resulting in fewer and delayed I-waves (Aberra et al., 2020; De Geeter, Crevecoeur, Leemans, & Dupré, 2015).

Di Lazzaro et al. (2012) has previously described a minimal cortical circuit hypothesised to generate I-waves. Components of this circuit were tested using the motor cortex model to establish the contributions to synchronised wave generation. A diagram of the proposed circuits by Di Lazzaro et al. (2012) are shown in Figure 6 below. To create these simpler circuits, the same neuron indices were used for the populations of layer 2/3E, 2/3I, and 5E and only the additional neuron groups and connections were removed so the network in the relevant populations were identical to the larger network model.

All differential equations were solved using a linear numerical integration with a step size of 0.1 ms. A total of 250 ms of neuronal activity was simulated, and the last 200 ms of steady-state model behaviour was analysed. The TMS activation was applied 150 ms into the simulation. The model took approximately 5 minutes to run on a six-core CPU. Firing frequencies were calculated from the instantaneous spiking activity of the entire population in each time window and smoothed using a Gaussian weighted average with a standard deviation of 0.15 ms.

The spiking frequency of neurons during each simulation was recorded, and I-waves were interpreted to be the cumulative activity of layer 5 excitatory neurons. The axons of layer 5 excitatory neurons made up the majority of the corticospinal tract where I-waves are physiologically recorded.

The model replicated the high-frequency activity observed in layer 5 pyramidal tract neurons in response to a simulated TMS. I-wave activity has been recorded from the epidural space of the cervical spinal cord in humans following simulation of the motor cortex. The measured changes in electrical potential were assumed to be the summation of the action potentials in the axons of layer 5 neurons, which extend down from the motor cortex (Esser et al., 2005; Rusu et al., 2014). From the simulation, a D-wave at the time of TMS delivery and series of three I-waves of neuron activity with interspike intervals (ISIs) of 1.2–1.5 ms was observed. These results were similar to the response recorded in humans and is shown in Figure 7A. A delayed two-wave response with lower amplitudes in the model was observed in response to simulated AP TMS as shown in Figure 7B.

Figure 8 shows different numbers of waves appear in each layer of the cortical model during PA stimulation. This model suggests that all cortical layers are involved in the generation of I-waves, but further analysis discussed later in this section suggests that particular groups may be key contributors. A raster plot of neuron firing arranged by firing frequency is shown in Figure 9. This plot shows that neurons contribute to the firing behaviour of the population unequally, a result supported by experimental measures of neuron firing in the cortex (Dabrowska et al., 2021; Lacey et al., 2014; Roxin, Brune, Hansel, Mongillo, & van Vreeswijk, 2011; Tomov, Pena, Zaks, & Roque, 2014). Also, just over 50% of neurons (6,566 out of 10,944) in the group contributed to the I-waves. Localised firing activity was observed during spontaneous activity, though during TMS a more widely distributed firing pattern was observed in the model (Figure 10).

A cortical silent period was observed immediately after an applied TMS, similar to those seen experimentally (Figure 11). The silent period is described as a pause in electromyographic activity following the onset of a motor-evoked potential in response to TMS during voluntary muscle contraction and is thought to be mediated by both cortical and spinal mechanisms, though the contributions of each are still unclear. Complete loss of the silent period has been observed in patients with ischemic lesions in the motor cortex, suggesting that the silent period is of cortical origin (Schnitzler & Benecke, 1994). Silent period duration shows large individual variation but typically ranges from 100–300 ms with duration increasing as TMS strength increases (Ahonen, Jehkonen, Dastidar, Molnár, & Häkkinen, 1998). Figure 11 shows a period of approximately 50 ms with a substantially decreased level of spiking activity following TMS, particularly in excitatory neurons, which supports the theory that cortical mechanisms may contribute to the silent period (Özyurt et al., 2019; Schnitzler & Benecke, 1994).

In the PA direction, the strength of TMS affected the amplitude of the resulting I-waves. Figure 12 shows the amplitude of waves in response to increasing the strength of PA TMS, showing a similar trend to that recorded by Nakamura, Kitagawa, Kawaguchi, and Tsuji (1996). Increasing the strength of the TMS, by increasing the proportion of neurons stimulated, increases the resulting amplitude of the population frequency response. The I3 wave was recruited at a higher TMS strength in this model, with the I3 wave only appearing when 22.5% of the neuron population was stimulated. This was similar to experimental results for the proportion of neurons affected by single pulse TMS to evoke an effect showing burst of action potentials in neurons recorded by Romero et al. (2019). The I3 wave is also thought to play a disproportionately significant role in the recruitment of motoneurons and size of the MEP (Thickbroom, 2011).

An increase in the strength of inhibitory connections (by 1.4 times) to simulate the effect of pharmacological GABA_A agonists resulted in an observed decrease in the amplitude of I2 and I3. Figure 13 shows the population-based firing activity following TMS when the strength of inhibitory connections was increased. This finding corresponds to observations in pharmacological experiments with TMS, where administration of a benzodiazepine such as lorazepam suppresses later I-waves but not the first I-wave (Di Lazzaro et al., 2012). Similarly, I1 waves are rarely affected by conditioning stimuli or enhancement of GABAergic activity (Di Lazzaro & Ziemann, 2013; Zewdie & Kirton, 2016).

Next, we assessed the relationships between the temporal parameters of individual neurons, refractory period, and transmission delay, and the amplitude and timing of I-waves. A decrease in refractory period led to a higher frequency response with the time between the waves decreasing from 2 ms at a refractory period of 1.5 ms, to 1.2 ms at a refractory period of 1 ms (Figure 14). Changes to excitatory and inhibitory delays had the biggest effect on later I-waves (i.e., I2 and I3), but in opposing directions. Increasing the excitatory delay resulted in a left shift (earlier in time) of the later I-waves, whereas increasing the inhibitory delay resulted in a right shift (later in time) of the later I-waves (Figure 15 and Figure 16). The amplitude of the D-wave remained constant, as expected, and I1 was also consistent except in some cases when the later I-waves may have superimposed, meaning activity of the I1 and I2 waves was combined, resulting in an increased amplitude.

Simplified circuits with specific components of the full model omitted were tested to investigate the contributions of different network components to the generation of I-waves. A minimal circuit that contained superficial layer 2/3 neurons and deep layer 5 neurons, as well as a group of inhibitory neurons based on previous models put forth by Di Lazzaro et al. (2012) were tested. The initial circuit with no recurrent connections between the layers (Figure 17) resulted only in an I1 wave being produced, suggesting that the direct connections between layer 2/3 and 5 are responsible for the generation of I1, and that more complex circuits with recurrent connections are necessary to produce later I waves.

I-waves have previously been hypothesised to be generated by recurrent cortical layer 2/3 neurons (Seo et al., 2016). However, the addition of recurrent connections to excitatory neuron groups without recurrent inhibition did not result in the generation of an I2 wave and only in combination with recurrent inhibitory connections were later I-waves produced (Figure 18). Inhibitory recurrent connections are necessary for an I2 wave to be produced, in addition to the initial I1 wave (Figure 19).

Recurrent inhibitory connections are vital for stable, balanced activity and the transient nature of the oscillations in response to stimulation. Increased excitation, through recurrent excitatory and inhibitory connections (which increase excitation and decrease inhibition in the neuron populations) such as in the circuit shown in Figure 20A, resulted in a higher rate of firing in the network and continuous synchronised oscillations (Figure 20B), in the model even in the absence of spontaneous input or stimulation.

Polysynaptic connections to and from other cortical layers such as 4 or 6 also contribute to the generation of later I-waves. Without layer 4 connections in the circuit, the model only generated I1 and I2; without layer 6 connections, the model generated three I-waves but with a lower amplitude I3 that, as mentioned previously, may be important for the recruitment of motor units and the generation of the MEP (Figure 21). These data suggest that while layers 2/3 and 5 are critical for I-wave generation, other cortical layers might also play a role in shaping the I-wave pattern.

The effect of paired pulse TMS on the motor-evoked potential, usually measured at a distal muscle by electromyography (EMG), can be either facilitatory or inhibitory of motor evoked potentials depending on stimulus strength and the interval period (Müller-Dahlhaus, Liu, & Ziemann, 2008). Short interval intracortical faciliation occurs when the first (conditioning) stimulus is, or both the conditioning and subsequent test stimulus are, administered at or above the motor threshold required to elicit motor-evoked potentials, and stimuli are administered around 1.5 ms, 3 ms, or 4.5 ms apart. Alignment of the interstimulus intervals with the timings of the I-wave peaks results in facilitatory interaction since the waves and stimuli are in-phase (Di Lazzaro et al., 2012; Di Lazzaro & Ziemann, 2013; Ziemann, Tergau, Wassermann, et al., 1998a).

Using our model of the motor cortex, paired pulse stimulation with an interstimulus interval of 2 ms, to match the interval of the I-waves, resulted in the increased amplitude of I1, I2, and I3 with an extra wave (i.e., a small I4) (Figure 22). Paired TMS at short intervals of 1.0 to 1.4 ms has been shown experimentally to be accompanied by larger and more numerous I-waves, with a more significant increase in I2 and I3, which matches our model results (Di Lazzaro et al., 1999; Ziemann & Rothwell, 2000). The observed increased I-wave amplitudes in the model could contribute to explaining short interval facilitation observed in paired pulse protocols, which increase the size of the resulting MEP when a second pulse is administered at specific time intervals due to the interaction of I-waves at the level of the motor cortex (Ziemann et al., 1998a).

With a homogonised, random connectivity in the model with no spatial connection topology but with preserved vertical connections between layers and the same number of synapses between neuron groups, a series of waves were still produced, though at different amplitudes to localised connections. Random connections resulted in a larger I2 and I3 wave (Figure 23). This suggests that the topology of the network is also a factor in I-wave generation.

A spiking neural network model of motor cortex circuitry has been developed and used to simulate and explain the cortical response of TMS on corticospinal neurons. The model replicates the input-output characteristics of single and paired pulse stimulation, as well as differences between PA and AP TMS, in regards to I-wave activity as synchronous high-frequency network activity. The effects of TMS strength and increased inhibition were also captured by this cortical circuit model. The model incorporated excitatory and inhibitory populations of individual neurons, replicating the known physiological structure of the cortex including layers and local connectivity. This circuit was also compared to simpler circuits previously proposed to be involved in the generation of I-waves (Di Lazzaro et al., 2012). Connectivity between neurons was defined at the population level, and consequently the resulting model gives insight into how groups of neurons can interact to produce cortical responses observed experimentally.

Di Lazzaro et al. (2012) hypothesised that the physiological generation of I-waves could be minimally represented by interactions between the circuits of inhibitory and excitatory neuron groups in layers 2/3 and 5 of the motor cortex. Our model uses simple LIF neurons situated in a large-scale cortical circuit based on previously published electrophysiological connectivity data to explain the origin of I-waves. In addition to the full model, we tested simpler circuits representative of those proposed by Di Lazzaro et al. (2012). The model demonstrated that I-waves can be explained by the recurrently connected circuitry of groups of neurons in the cortex. The results of this model support the idea that, while a subset of the connectivity is responsible for the generation of the initial I-waves, other cortical layers are likely involved in the generation of the full pattern of produced I-waves, rather than simply layer 2/3, and 5 which have been the main focus of previous work (Di Lazzaro et al., 2012; Rusu et al., 2014).

The construction of this cortical model contained excitatory and inhibitory neuron groups in layers 2/3, 4, 5, and 6 with specific connections within and between each layer based on neuroanatomical data (Binzegger et al., 2004; Thomson, West, Wang, & Bannister, 2002). An important advancement of our model beyond that of Esser et al. (2005) is the inclusion of a layer 4 neuron group, which has recently been shown to exist in the motor cortex (Barbas & Garcia-Cabezas, 2015; Yamawaki, Borges, Suter, Harris, & Shepherd, 2014). Interestingly, the model reproduced three I-waves without the layer 6 excitatory and inhibitory neuron groups, but not without layer 4 neuron groups, which differs to the Esser et al. (2005) model, which represented layers 2/3, 5, and 6. The difference between the layers’ effect on the generation of I-waves may be due to the differing density of connections between the groups with more dense connections in layer 4. We took a similar modelling approach to Esser et al. (2005) but used different neuron models and connectivity, though both connectivities were physiologically based. Both spiking neural network approaches were able to reproduce I-waves, suggesting that recurrent, polysynaptic circuitry in the cortex is critical for generating synchronised oscillatory activity.

Paired pulse protocols can result in either the facilitation or the inhibition of the motor evoked potential Di Lazzaro et al. (1998). Epidural spinal cord recordings have shown larger and more numerous I-waves with paired-pulse TMS at short interstimulus intervals of 1.0–1.4 ms, demonstrating the intracortical origin of short-interval cortical facilitation (SICF), as replicated in this model (Di Lazzaro et al., 1999). A subthreshold stimulus, delivered more than 6–25 ms prior to a suprathreshold stimulus, resulted in facilitation whereas a subthreshold stimulus delivered less than 5 ms followed by a suprathreshold stimulus suppresses the MEP amplitude (Peurala, Müller-Dahlhaus, Arai, & Ziemann, 2008; Zewdie & Kirton, 2016). Suppression of late I-waves has also been observed by direct epidural recordings, suggesting a subthreshold stimulus may activate intracortical inhibitory circuits (Di Lazzaro et al., 1998). A subthreshold paired-pulse stimulus activating inhibitory neurons could be explored further, with changes to the TMS and inhibitory neuron model. In addition, an inhibitory-specific neuron model and the resulting changes to network dynamics as well as I-wave generation could be future areas of investigation.

GABAergic interneurons play a key role in modulating excitatory activity within circuits of the cortex. GABA agonists have an effect on the amplitude and increased inhibition can suppress later I-waves, as demonstrated in this model (Tremblay, Lee, & Rudy, 2016; Ziemann, Tergau, Wischer, Hildebrandt, & Paulus, 1998b). This model also illustrated the importance of recurrent inhibitory connections in I-wave generation, particularly in the generation of I2. Inhibitory interneurons are thought to have important roles in cortical functions, including modulation of cortical circuits through gain, regulation of pyramidal cell firing, plasticity, synchronisation, and the generation of cortical patterns necessary for information transfer (Tremblay et al., 2016). This model shows that inhibitory connections are important in the generation of I-waves and corticomuscular transmission.

The differing results generated by testing variations in circuitry in this cortical model supports the idea of the structure of neural circuitry contributing to functional activity. Tsubo, Isomura, and Fukai (2013) suggested that the layers work as functional units with layer 2/3, serving as a “resonant oscillator” and layer 5 pyramidal neurons operating as “integrators.” Other hypothesised mechanisms of I-wave generation include the individual pyramidal cells of layer 5 operating as neural oscillators, repetitively discharging in response to depolarisation with high-frequency bursting due to their intrinsic membrane and geometric properties (de Kock et al., 2021; Hedrick & Waters, 2011; Kernell & Chien-Pingt, 1967; Silva, Amitai, & Connors, 1991; Triesch, Zrenner, & Ziemann, 2015). Neurons anatomically span several layers and are typically organised into layers based on the location of their cell body (Larkum, Petro, Sachdev, & Muckli, 2018). Although the laminar organisation in the cortex was preserved in this model, the heterogeneity of neurons in regards to morphology and electrophysiological behaviour even within each layer cannot be understated (Oswald, Tantirigama, Sonntag, Hughes, & Empson, 2013). Specifically in the motor cortex, layer 5 has previously been separated into layer 5a and 5b due to the cells being larger and more numerous in the lower section of layer 5 (Sherwood & Hof, 2007).

Our representation of homogeneous populations of single-compartment neurons in excitatory and inhibitory neuron groups cannot assess the consequence of neuron morphology, which plays a role in the integration of input at the neuron level (Seo et al., 2016). More detailed compartment models of layer 5 pyramidal neurons have been previously developed by Rusu et al. (2014), whereas our work aims to describe possible roles of connectivity and neural population dynamics in the generation of I-waves and, for practical computational reasons, benefited from modelling simple point neurons. The complexity of simulating many thousands of multiunit neuron models is still limited by computational power. Previous models have also explored complex neurotransmitter and synaptic models (Esser et al., 2005; Rusu et al., 2014), whereas our model used a simpler neuron model, each with only a single type of excitatory and inhibitory synapse. A recent review by Ziemann (2020) stated that there is currently no evidence that I-waves are generated by high-frequency oscillations of individual corticomotoneuronal cells, favouring models of polysynaptic excitatory and inhibitory interneuronal circuits. Although more detailed models could be introduced into our modelling framework, the ability of a simpler population-based model to replicate experimental results suggests that complex neuron-level behaviour may not be essential for the generation of I-waves.

Recent work in intracortical recordings of the motor cortex suggests a population-based, dynamical systems approach to understanding neural activity (Kaufman, Churchland, Ryu, & Shenoy, 2014; Shenoy, Sahani, & Churchland, 2013; Yu et al., 2005). Neurons of the motor cortex neurons are shown to respond even in the absence of movement, for example, during motor imagery or preparatory activity. Therefore, it is unlikely that activity of individual neurons are able to provide functional links to an organisms’ behaviour (Churchland, Cunningham, Kaufman, Ryu, & Shenoy, 2010; Facchini, Muellbacher, Battaglia, Boroojerdi, & Hallett, 2002). This model, which is on the scale of representing neural populations in the cortex, can begin to link the structures of cortical circuits to functional output. Input structures such as premotor cortex and basal ganglia to drive cortical activity will also be necessary (Logiaco, Abbott, & Escola, 2021). To further investigate the functionality of the cortex using this model, it could be placed within a larger scale framework of motor generation to connect input, motor cortical activity and resulting muscle force generation (Haggie et al., 2023b).

The model used a simplified version of TMS, stimulating a random selection of neurons making up a defined proportion of the population. However, the sites of real TMS stimulation are shown to be dependent on individual neuron properties, such as orientation of the axons relative to the electric field (Pashut et al., 2014; Rotem & Moses, 2006, 2008; Siebner et al., 2022). Previous models of TMS have included computed electrical fields and volume conductor models, which take into account complex geometry of the cortex including gyral folding (Seo et al., 2016). Computational models show that the geometry of cortical folding as well as neuronal morphology are key parameters of the effects of TMS (Seo et al., 2016). For the purpose of explaining the origin of I-waves, detailed geometry for stimulation was not taken into account and more complex modelling of the TMS effect on the cortex, for example, as an electromagnetic field stimulating axons lying in specific orientations, rather than as direct current inputs to point neurons, could be a future area of exploration.

Future work building on this research could extend the model to include spinal and neuromuscular components to produce motor evoked potentials, further exploring the cortical response and the downstream effects of various experimental TMS protocols. Next steps could also explore the effect of paired pulse stimulation on inhibition, which could suppresses late I-waves in short-interval cortical inhibition experiments. It is suggested that subthreshold TMS mostly activates inhibitory neurons because of their lower threshold for magnetic stimulation (Rusu et al., 2014). Higher stimulation levels might also activate more excitatory neurons or possibly neurons from both AP and PA circuits, contributing a higher number of I-waves being generated as observed in human experimental data (Aberra et al., 2020; Siebner et al., 2022).

Using a spiking neural network model of motor cortical circuits, we have demonstrated how populations of neurons can reproduce observed responses to TMS. The model recreates observations in spontaneous firing, I-wave activity, and responses to increased inhibitory activity in the motor cortex. We also demonstrated how different circuits could contribute to the effects of AP and PA direction TMS, and tested previously hypothesised circuits thought to be involved in the generation of I-waves to show the necessary contribution of more complex connectivity, which exists physiologically in the cortex. Spiking neural network models involving large numbers of individual cells are advantageous over morphological models or mean-field models because they enable the observation of contributions of both individual neurons and population-based activity, linking between micro- and mesoscale physiological activity. Developing computational models alongside experimental testing is vital for contributing to the understanding of the effects of TMS and the contribution of cortical circuits to the generation of movement.

The cortical response to TMS can be explained by recurrent connections in a spiking neural network model of cortical circuitry not tuned to generate oscillatory responses. Parameters that played a significant role in shaping I-wave generation include the refractory period and the transmission delays of neurons. The model explained the generation of I-waves in populations of neurons of the cortex by implementing and testing a range of previously defined circuits. Our results verify that layer 2/3 to layer 5 connections are important in generating I1, and recurrent inhibitory connections play a role in generating I2. More complex circuitry involving other cortical layers was necessary to achieve the generation of I3 as well as realistic resting-state activity. With this model, we also proposed a circuit to reproduce the AP response and showed the effects of specific model parameters and increased inhibitory activity. This work extends previous studies of cortical modelling to understand experimentally measured behaviour in the human nervous system.

The author thanks John Cirillo for their valuable discussions and Mark Sagar for initial project funding and administration.

Lysea Haggie: Investigation; Methodology; Visualization; Writing—Original draft. Thor Besier: Supervision; Writing—Review & editing. Angus McMorland: Conceptualization; Supervision; Writing—Review & editing.

Transcranial magnetic stimulation:

A noninvasive technique for brain stimulation involving a fluctuating magnetic field which induces electrical currents in the neurons of the cortex.

I-waves:

A repetitive oscillatory response measured in the spinal cord which occurs in response to transcranial magnetic stimulation over the motor cortex.

Spiking neural network:

An artificial neural network that mimics properties of biological neurons such as the membrane potential, threshold, and action potential dynamics.

GABA:

The main inhibitory neurotransmitter in the central nervous system that reduces the neuron membrane potential away from firing threshold.