Numerical Modeling of Plasticity Induced by Quadri-Pulse Stimulation

Quadri-pulse stimulation (QPS), a type of repetitive transcranial magnetic stimulation (rTMS), can induce a considerable aftereffect on cortical synapses. Human experiments have shown that the type of effect on synaptic efficiency (in terms of potentiation or depression) depends on the time interval between pulses. The maturation of biophysically-based models, which describe the physiological properties of plasticity mathematically, offers a beneficial framework to explore induced plasticity for new stimulation protocols. To model the QPS paradigm, a phenomenological model based on the knowledge of spike timing-dependent plasticity (STDP) mechanisms of synaptic plasticity was utilized where the cortex builds upon the platform of neuronal population modeling. Induced cortical plasticity was modeled for both conventional monophasic pulses and unidirectional pulses generated by the cTMS device, in a total of 117 different scenarios. For the conventional monophasic stimuli, the results of the predictive model broadly follow what is typically seen in human experiments. Unidirectional pulses can produce a similar range of plasticity. Additionally, changing the pulse width had a considerable effect on the plasticity (approximately 20% increase). As the width of the positive phase increases, the size of the potentiation will also increase. The proposed model can generate predictions to guide future plasticity experiments. Estimating the plasticity and optimizing the rTMS protocols might effectively improve the safety implications of TMS experiments by reducing the number of delivered pulses to participants. Finding the optimal stimulation protocol with the maximum potentiation/depression can lead to the design of a new TMS pulse generator device with targeted hardware and control algorithms.


I. INTRODUCTION
Transcranial magnetic stimulation (TMS) is a non-invasive method for modulating central nervous system functions. Brief stimuli of an electromagnetic field are applied to the neural tissues with current-carrying treatment coils. The fast-changing field induces electrical currents within the neurons and modulates cell activity. Certain paradigms for magnetic stimulation can increase or decrease cortical excitability even after the end of stimulation which is called plasticity. The induced plasticity generally depends on the time interval between pulses, the stimulation intensity, the session duration The associate editor coordinating the review of this manuscript and approving it for publication was Taous Meriem Laleg-Kirati . and the waveform shape (monophasic, biphasic) [1]. Certain patterns can strengthen or weaken synapses between neurons, a process called long-term potentiation (LTP) and depression (LTD), respectively. These effects generally outlast the period of magnetic stimulation and the aftereffects are highly dependent on the stimulation characteristics.
The understanding of how magnetic stimuli interact with cortical circuits is improved by the application of biophysically-based models which commonly characterize neurophysiological mechanisms with equations. For instance, field-based modeling of the brain activity is now well established in neural field or neuronal population models, where neural dynamics are modeled in terms of axonal flux rates of neural populations and mean firing rates, rather than modeling individual neuron dynamics [2]- [4]. Population models are particularly suited to the modeling of TMS since a pulse can stimulate several cm 2 of cortex. Various neural populations, such as excitatory neurons or inhibitory interneurons, can be connected together, with predefined parameters and weights as the coupling strength between populations.
Plasticity can be added into the neural field modeling by spike timing-dependent plasticity, in short STDP. This model describes biophysical processes with equations for the firing rates of excitatory and inhibitory populations of stimulated neurons, where the time interval between presynaptic and postsynaptic spikes establishes the alteration in synaptic strength. The key parameters of this model are neural firing rates, axonal pulse rates and the mean soma potentials [5]. Broadly, the average firing rate of cells is a function of the average spike rate in axons and the average membrane potential. Another approach to incorporate the plasticity into neuronal population-based modeling is Calcium-dependent plasticity (CaDP) theory [6]. The CaDP represents cellular calcium dynamics and its consequences on the synaptic strength. Typically, moderate concentrations of calcium lead to depression, while high concentrations lead to potentiation [7].
A new TMS protocol that has not yet been reproduced by the neural plasticity models is the quadri-pulse stimulation (QPS). Despite the early promise of repetitive biphasic-pulse protocols, such as theta-burst stimulation (TBS), later research reported considerable variability between individuals, where the same rTMS modality might induce opposite plasticity effects for different participants [8]. One novel paradigm that is showing promising aftereffects is the QPS, which uses monophasic stimuli instead [9], [10]. This method utilizes four monophasic stimuli at predefined inter-stimulus intervals (ISIs) in a burst with an inter-burst interval of 5 s, where modulation of cortical excitability highly depends on the ISIs [11]. According to human experiments, recorded motor-evoked potentials (MEPs) are weakened for QPS50 (ISI = 50 ms), whereas they are strengthed for several hours after QPS5 (ISI = 5 ms) [12]. The subsequent experiments also showed that short ISIs have (LTP)like aftereffects (i.e. QPS1.5, QPS5, and QPS10), while long ISIs have (LTD)-like aftereffects (i.e. QPS30, QPS50, and QPS100) [11], [10].
Producing short ISIs is one of the main challenges of TMS devices. For example, in the majority of literature, to generate the QPS protocol, four Magstim TMS 200 devices (Magstim, UK) are connected with a specially designed combining module [11]. Therefore, it increases the size and cost of the required QPS system. To address this need, a new TMS device, called controllable TMS or cTMS, was introduced by Peterchev et al. [13]. Using new power electronic switches and circuit architectures, the cTMS is able to generate consecutive stimuli with a minimum ISI of 2 ms, as well as change the pulse width (i.e. the pulse waveform). For the QPS protocol, the cTMS can generate unidirectional pulses, instead of conventional monophasic pulses. The induced electric field on the cortex by unidirectional pulses has three different phases: a negative, a positive and again a negative phase. The dominant part of this stimulus is a positive phase and it acts almost like a monophasic pulse [14].
The new studies emphasize that there are still many stimulation parameters that can be controlled that have not yet been fully investigated, which might become a future basis to optimize or invent novel magnetic stimulation paradigms. Using the concept of pulse width modulation, Memarian et al. were able to synthesize arbitrary waveforms [15]. These designs open a novel parameter space for TMS by generating and manipulating a magnetic stimulus whose stimulus shape has been effectively controlled. The pulse width and direction change can help to selective target neuronal populations (so-called neural tuning) [16].
Furthermore, one of the unsolved challenges in plasticity modeling is the effect of the magnetic pulse shape on the TMS-induced plasticity, because, in conventional TMS devices, the waveform alteration was limited. According to recent studies, to introduce and develop a new TMS device for the generation of flexible pulses and new rTMS protocols with different waveforms [17]- [19], the importance of further research in the field of neural plasticity modeling is growing. The cTMS device is one of the most important neurostimulators in the TMS field which enables the test of different stimulus width to reach the highest plasticity.
Guided by existing plasticity models and experimental QPS results, we modeled QPS-induced plasticity based on neural population models and the STDP model, as illustrated in Figure 1. Additionally, the possibility of implementing QPS protocols with unidirectional pulses and the effects of different pulse widths have been explored for the first time.

II. METHOD
To find changes in synaptic weights between two populations of neurons, we use the phenomenological STDP approach, where the synaptic weight alteration between populations is dependent on the time difference between the post-synaptic and pre-synaptic spikes. We base this study upon the well-documented STDP implementation of Wilson et al. in [20] and the foundation provided by [1], [5] to modeling neuronal populations and cortical activity to capture the plasticity response. The key principles are as follows. The average characteristics of a population of neurons are assumed, instead of modeling the behavior of a single neuron. The cortex is modeled by a group of excitatory neurons coupled with a group of inhibitory cells. To include the effect of the TMS pulse, a third population of excitatory cells was added which connects to both excitatory and inhibitory neural populations. Figure 1 displays the described model. The main parameters of STDP description include the following mathematical modeling: alternation of the target synaptic coupling, mean firing-rate of the axonal input, axonal spike rate, frequency response of the excitatory cells, axonal transfer function, and frequency spectrum for the magnetic stimuli.

A. QUANTIFYING THE NEURAL PLASTICITY
The mean relative change of the target synaptic coupling ( W ee ) (in the case of excitatory-to excitatory synapses) was calculated within 10 seconds of a particular TMS protocol and normalized by the number of stimuli. This is expressed as: where a TMS protocol repeats over a time-scale T (here T = 10 s), N T is the total pulse number, W ee (t = 0) is the value of the synaptic weight at the starting point of the protocol or initial plasticity (baseline) and the angle brackets depict a temporal average over t. This equation enables the identification and classification of protocols in terms of potentiation (increasing in strength or W ee > 0) or depression (decreasing in strength or W ee < 0). Synaptic weight changes over time (dw ee /dt) are calculated by the following equation [20]: where Q e is the mean firing-rate of the axonal input, as a function of the average axonal spike rate (ϕ), and |Q e (ω)| 2 is the frequency response of the excitatory cells. e (ω) is the axonal transfer function and explains the reaction of the axonal firing rate to alterations in the firing rate Q e . H ee (ω) describes the STDP function, which depends on the post-synaptic neuron.
Re[ e (ω)H ee (ω)] is called the plasticity function. H ee (ω) is calculated by [5]: where A + e is the STDP positive weight constant and A − e is the STDP negative weight constant which are both dimensionless. t + e and t − e are the STDP positive decay constant and the STDP negative decay constant, respectively.
Mean firing-rate or Q e (ω) is calculated according to the following equation [5]: where G e and G i are the excitatory and inhibitory gains of the system respectively, where the gain, which is dimensionless, is the product of the mean synaptic weight, the average number of connections between cells and the rate of change of mean firing rate with respect to voltage. ρ TMS (ω) is the frequency spectrum for the QPS trains, λ ie and λ ee depict the strengths of the couplings between the ρ TMS and the inhibitory and excitatory groups, respectively. The transfer functions and variables used in the numerical modeling are shown Figure 1. L(ω) denotes the impulse response function of the soma potential of neurons responding to the incoming axonal input. For excitatory synapses: where the rate constants α ee and β ee describe the rise and decay rates of the soma response to the axonal input. These variables characterize the combined rates of the synaptic and soma dynamics. In the case of inhibitory synapses: This equation includes the effects of both gammaaminobutyric acid (GABA A and GABA B ) neurotransmitters, which are the chief inhibitory compounds in the central nervous system. In (6), the A and B subscripts describe the response to GABA A and GABA B respectively. α A is the neurotransmitter growth rate in response to the inhibitory GABA A input, α B is the growth rate in response to the inhibitory GABA B input, β A is the neuron decay rate in response to the inhibitory GABA A input and β B is the neuron decay rate in response to the inhibitory GABA B input. The axonal transfer functions for excitatory and inhibitory synapses are: 26486 VOLUME 9, 2021 FIGURE 2. Plots of the estimated and measured W ee percentage for a conventional monophasic pulse. Values greater than 100% denote long-term potentiation, values below 100% long-term depression. The red plots show the relationship between the measured mean MEP amplitude in human experiments and the inter-pulse interval. The red shading depicts the standard uncertainty in the mean, according to experimental results presented in [11].
where γ e and γ i are the excitatory axonal rate constant and the inhibitory axonal rate constant, respectively. The values for the standard parameters have been chosen from [20]. They are either measured experimentally or are fitted to the related physiological phenomena.

B. SOLVER
The described model is solved with algebraic symbols in Fourier space. The total time of the assumed QPS protocol is 10 seconds, which is discretized with a sampling frequency of 100 MHz. The total number of data points for each QPS protocol is equal to one million samples.

A. VALIDATION OF MODEL
To demonstrate the ability of the STDP and the neural field model to predict the aftereffect of the QPS on the plasticity, the described model was run with the conventional monophasic QPS protocols [11]. As far as we have studied, that research is the most comprehensive database related to the QPS protocol. By examining the stimulus-response function of QPS-induced plasticity in [12], [21], it can be seen that the results of their human trials almost overlap with the outcomes reported in [11]. The electrical circuits of the Magstim 200 neurostimulator were simulated in MATLAB Simulink software (Powergui blockset). A total of six QPS protocols with different time intervals between pulses (QPS1000-QPS1.5) were tested. The predicted results of the plasticity and the results of the human experiments, along with the monophasic pulse shape generated in the Magstim 200 device are shown in Figure 2 and Table 1. The results of the predictive model broadly follow what is typically seen in human experiments, although the timings of the LTD and LTP regions are not reproduced. precise timings of the LTD and LTP regions are not reproduced.

B. ESTIMATION OF PLASTICITY FOR UNIDIRECTIONAL PULSES
To estimate the effect of different waveforms on the neural plasticity, the near-rectangular pulses generated by the cTMS device have been explored. For this purpose, four unidirectional waveforms with positive phase widths of 120 µs, 100 µs, 60 µs and 40 µs, with different negative phases (16 waveforms in total) have been generated by circuit simulation. These results and waveforms are shown in Figure 3 (a)-(d). Specifically, stimulation at ISIs of more than 31 ms (QPS30, QPS50 and QPS100) induced potentiation, and stimulation at ISIs of less than 31 ms (QPS2, QPS5 and QPS 10) induced depression. Due to the limitations of the cTMS device, it was not possible to generate the QPS1.5 protocol, so the nearest feasible protocol (QPS2) was explored.
By comparing the plasticity results for conventional and near-rectangular pulses, it is found that unidirectional pulses can induce almost the same plasticity. However, the highest difference is seen in the QPS10 protocol. Also, as the width of the positive phase increases, the value of potentiation will increase. This increase is evident in the QPS5. Furthermore, in pulses with the same positive width, each pulse with a smaller negative width induces more potentiation. According to the simulations, the effect of pulse waveforms on depression was limited.

IV. DISCUSSION
In this study, we present how models of neural populations and physiologically-based theories of plasticity can imitate experimental results and be used to test predictions for new TMS pulses. For this purpose, we implemented the QPS-induced plasticity scheme together with the two-population neural field model and the STDP model in [5], [20]. We studied its theoretical properties and compared the model results with existing experimental results.
In particular, we show that the previously neglected magnetic pulse shape plays a role in plasticity modeling, in a total of 117 different scenarios. By combining the TMS circuit and the plasticity models, the optimal stimulation protocol for the maximum plasticity can be determined.
For the first time, the proposed algorithm has been utilized to estimate the plasticity induced by the QPS protocol. The key innovation in this research is the use of the frequency spectrum of the applied TMS pulses in the model while in numerical models in literature the stimulus is assumed to be in the form of an ideal pulse. For instance, the stimulus waveform is modeled as a delta-function impulse in [20], a top-hat function of length 0.5 ms in [6] and a 500 µs rectangular pulse in [22]. Evidently, those models can not investigate the effect of new waveforms on the plasticity.
The results of QPS-modeling reveal that firstly, unidirectional pulses can achieve a similar effect in terms of TMSinduced plasticity, secondly, pulses with wider positive phase and narrower negative phase can induce stronger plasticity. It is predicted that if unidirectional pulses are utilized in the QPS protocol, the most potentiation will be seen in the QPS5 protocol and the most depression will be observed in the QPS50 protocol, as shown in Fig. 3.
By comparing the error values of the validation step in [20], related to the paired-pulsed modeling results and this work, it can be concluded that using the real spectrum of magnetic stimuli instead of modeling it with ideal pulses may reduce the computational error and bring the simulation closer to the human experiments.

A. LIMITATIONS
The main restriction in implementing the proposed model is to find the accurate parameters for the dynamics of a neuron and the plasticity-related equations (3)(4)(5)(6)(7)(8), since the equations' parameters might nonlinearly fluctuate with gender, age and race. Our study, like all biophysically-based research, is sensitive to the choice of parameter values. Indeed, biophysical variation in parameter values may account for the wide variation in experimental results for some forms of TMS. In this work, we used the standard parameter values from the plasticity modelling of [20] without further optimization to the QPS results. Notable matching of our results with the outcomes of the human QPS trials shows that the parameter choice is reasonable and that any sensitivities are not likely to be greatly significant. Plausible parameter values and sensitivity of neural field modelling to them have been discussed in detail in [23].
Additionally, the neural field modeling is extracted from physiological principles at a microscopic level, while the TMS device operates at a macroscopic level, both in terms of applying the stimulation and the evaluation of the results, using electroencephalography or MEP measurements [24].
In summary, finding realistic parameters for modeling of plasticity neural behavior, and TMS-modeling at the microscopic level, instead of macroscopic level, are key error sources between experimental results and modeling outcomes, especially the differences seen in the QPS10 protocol. Therefore, in future studies, focusing on these limitations and a mapping between microscopic-macroscopic levels should be considered to achieve better accuracy.
On the other hand, the obvious differences between the frequency spectrum of the unidirectional pulses, and the conventional monophasic pulses, can be another factor in the differences observed in the induced plasticities, which may emphasize the effect of magnetic pulse waveform on the neural behavior.

V. CONCLUSION
Quadripulse stimulation has recently been introduced to induce more permanent plasticity. We explored a biophysically-informed model based on the spike timingdependent plasticity and compared this with the results of existing experiments. The plasticity estimation allows researchers to optimize the stimulus waveform for controllable magnetic pulses, especially for the cTMS device. This prediction can reduce the number of tests needed to find the optimal plasticity for unidirectional cTMS pulses. Thus, the introduced algorithm demonstrates a useful mechanism to optimize TMS-induced plasticity.