key: cord-0799289-59wc8gln
authors: Memarian Sorkhabi, Majid; Wendt, Karen; Wilson, Marcus T.; Denison, Timothy
title: Estimation of the Motor Threshold for Near-Rectangular Stimuli Using the Hodgkin–Huxley Model
date: 2021-05-31
journal: Comput Intell Neurosci
DOI: 10.1155/2021/4716161
sha: aeb62f7fd4204b85c81ef8fa99eabfb1ec7166b0
doc_id: 799289
cord_uid: 59wc8gln

The motor threshold measurement is a standard in preintervention probing in TMS experiments. We aim to predict the motor threshold for near-rectangular stimuli to efficiently determine the motor threshold size before any experiments take place. Estimating the behavior of large-scale networks requires dynamically accurate and efficient modeling. We utilized a Hodgkin–Huxley (HH) type model to evaluate motor threshold values and computationally validated its function with known true threshold data from 50 participants trials from state-of-the-art published datasets. For monophasic, bidirectional, and unidirectional rectangular stimuli in posterior-anterior or anterior-posterior directions as generated by the cTMS device, computational modeling of the HH model captured the experimentally measured population-averaged motor threshold values at high precision (maximum error ≤ 8%). The convergence of our biophysically based modeling study with experimental data in humans reveals that the effect of the stimulus shape is strongly correlated with the activation kinetics of the voltage-gated ion channels. The proposed method can reliably predict motor threshold size using the conductance-based neuronal models and could therefore be embedded in new generation neurostimulators. Advancements in neural modeling will make it possible to enhance treatment procedures by reducing the number of delivered magnetic stimuli to participants.

Transcranial magnetic stimulation (TMS) is a noninvasive and promising tool to modulate the human cortex and can be used to induce lasting changes in neural activity, with both research and clinical applications [1, 2] . FDA has cleared TMS for the treatment of various neurological and psychiatric diseases and it is under consideration for many other disorders. In TMS devices, strong, brief current stimuli driven through a wire-wound coil induce an electric field to modulate nerve cells. Despite the success of the TMS procedure, there are several critical limitations associated with the magnetic pulse shape parameters. Conventional neurostimulators most commonly employ LC resonant circuits, resulting in a constant waveform for the magnetic pulses, depending on the predetermined hardware [3, 4] .

More adjustable control of the magnetic field waveform can potentially enable new research and clinical applications that are not realizable with conventional TMS tools, such as changing the waveform of the magnetic pulse [4] . Addressing this need, Peterchev et al. have introduced a controllable-TMS (cTMS) machine to produce flexible near-rectangular stimulus shapes with three different architectural variations [5] . Utilizing the concept of pulse width modulation (PWM), Sorkhabi et al. can synthesize arbitrary waveforms [6, 7] .

ese Instruments open a new parameter space for magnetic stimulation by producing and manipulating a magnetic pulse whose pulse shape has been impressively managed. With the control of PW and pulse direction, neural populations could be selectively targeted [8] .

e motor threshold (MT) is one of the standard quantifiable measures of corticospinal excitability in TMS research and is used to determine the TMS intensity. e MT is usually defined as the lowest stimulation intensity able to produce motor-evoked potentials (MEP) of size 0.05 mV with a 50% probability in the abductor pollicis brevis [9] . e MT determination might be measured in the relaxed muscle (resting MT or RMT) or during active muscular contraction (AMT). Estimating the MT can beneficially improve the safety implications of preintervention probing in the magnetic stimulation tests by minimizing the number of delivered pulses to individuals [10] .

is modeling also allows the optimization of stimulation waveforms in terms of the required minimum hardware/energy to generate them. As a result, the possibility of optimizing simulation-based tools will be enabled.

In the current work, we simulate the cTMS circuit to generate different near-rectangular pulses, directions, and intensities. Four insulated-gate bipolar transistors (IGBTs) switches and freewheeling diodes, which form the two half-bridges structure, were used to connect the treatment coil to the DC pulse capacitors. Two isolated DC supplies (V DC1 , −V DC2 ) have been used to charge the DC pulse capacitors and the output stimulus can have four different voltage levels (

e cTMS circuit structure with the parameters presented in [11] was simulated in the MAT-LAB Simulink environment (Powergui blockset, R2020a). e circuit model of this device and its typical parameters are shown in Figure 1S , supplementary data. e simulated results were used in the Hodgkin-Huxley (HH) model to estimate the MT.

e HH equations have been widely used in the modeling of several neuromodulation modalities and neuronal behavior [12] [13] [14] . In this study, for the computational modeling of cortical neurons, a Hodgkin-Huxley type model is selected, including four conductances (the voltage-dependent sodium, potassium, slow potassium, and resting (leak) membrane conductance) and the TMS-induced current density.

is model has biophysically meaningful and measurable parameters and allows researchers to investigate the interplay between an external stimulus current and an ionic current. To establish this conductance-based neuronal model, we review experimental data from various preparations and obtain a computational model that captures the fundamental aspect of the intrinsic characteristics utilizing a minimal number of equations. e HH model is more realistic than the RC integrated-fire model and can be fit for physiological measurements. e main equations and parameters are derived from [15] [16] [17] and the required parameters to apply the magnetic stimulation are added according to the following equations:

e C m term represents the membrane lipid bilayer capacitance. e symbols J Na , J K , J M and J leak denote the current densities of sodium, potassium, slow potassium (delayed rectifier and slow nonactivating current), and other leakage ions, respectively, whereas J stim is the TMS-induced current density by the cTMS stimulus (units: µA/cm 2 ). g Na , g K and g M are the maximal values of the conductances for the voltage-dependent sodium, potassium, and slow potassium currents per unit area. e dimensionless gating variables n, m, and h are associated with potassium channel activation, sodium channel activation, and sodium channel inactivation, respectively, and are bounded between 0 and 1. e transition rates between open and closed states of the ion channels (α i and β i , i � m, n, h) are positive coefficients that nonlinearly depend on the membrane potential V m , as determined below:

2

Computational Intelligence and Neuroscience

where p ∞ and τ p are the steady-state activation and time constant for the slow potassium current, respectively. As seen in equations (3)-(8), the conductances of the sodium and the potassium channels (g Na , g K and g M ) are a nonlinear function of time and the membrane voltage V m , while this parameter is constant for the leak channels. e details of the variables and their values involved in these equations are provided in Table 1 .

e equations of the single-compartment Hodgkin-Huxley type model were solved in MATLAB with Euler's numerical solution method (timestep: dt � 1 μs, number of iterations: 20000). To find a relationship between the treatment coil voltage (stimulation intensity) and the corresponding induced current in the cortex, the finite-element method (FEM) and 3D MRI-derived head model are used to compute the induced Electric-field (E-field). e head tissue of a healthy subject is assumed to be homogeneous with an electrical conductivity of 0.27 S/m [18] . E-field distributions were calculated with COMSOL Multiphysics 5.3 for a Magstim 70 mm figure-of-8 coil (P/N 9925, 3190). e specifications of the figure-8 coil are derived from [19] , which has 9 windings in each lobe with inner and outer diameters of 5.2 and 8.8 cm, respectively. e final mesh structure of the head comprised 120,000 nodes and 2.4 million tetrahedral elements. Figure 1 shows E-field norm distributions on planes 5, 10, and 15 mm beneath the coil, where the coil-to-scalp distance is 2 mm and the peak coil voltage is 1 kV. Given that the target hand motor area for finding the MT is generally assumed to be 15 mm below the scalp, it can be concluded from the FEM results that for every thousand volts of a figure-8 coil voltage, an electric field of 100 V/m is induced in the motor hand knob [20] .

erefore, the induced current to the target area will be 27 A/ m 2 per kV (J stim in equation 1). e results of human tests that have previously shown the relationship between the MT and the conventional monophasic and biphasic stimuli [21] [22] [23] have been used to validate the HH equations and ratios obtained by the FEM modeling. It was found that the FEM method and reported experimental results are consistent with an accuracy of ±4.6%. e approximation of the neuron's membrane as a low-pass RC filter and experimentally measuring the membrane potential changes with a pickup coil also shows similar results. For more details of this measurement, see [6] .

We set the transmembrane current density driving the HH model to this value. While transmembrane current densities are in practice not equal to the induced current densities, being dependent on cellular and brain geometry and conductivity variations, we expect them to be strongly related. We identify motor threshold as being the pulseamplitude required to produce an action potential in the HH model for a given pulse shape. Resting motor thresholds are consistently higher than AMTs, usually by about 10% of the stimulator output [24] . ere is also a positive correlation between participant's resting and active MTs [25, 26] .

erefore, by finding one of these two, the other value can be estimated.

To validate the estimated MT, we used the results of human experiments reported by [5, 8, [27] [28] [29] , in which 50 individuals were exposed to the different magnetic pulses of the cTMS neurostimulator. e data was collected from three types of cTMS waveforms, monophasic, unidirectional, and bidirectional in posterior-anterior (PA) or anterior-posterior (AP) directions. In those studies, the relationship between the pulse width and the effectiveness of a magnetic pulse (i.e., the strength-duration relationship) is specifically considered. A summary table indicating the number, population age, gender, and health conditions of subjects are available in Table 1S , supplementary data.

In [5, 27] , the effect of three monophasic stimuli (pulse width, i.e., positive phase duration, of 30, 60, 120 μs) on the RMT and the AMT was examined. ese pulses are shown in Figure 2 (a) and the measured RMT and AMT are displayed in Figures 3(a) and 3(b), respectively. As reported in those studies, by increasing the PW, the intensity of the stimulus required to see the defined MEP decreases. Also, the RMT is approximately 10% higher than the AMT.

In [8] , the effect of three unidirectional stimuli (negative pulse width of 40, 80, 120 μs) on the RMT and the plastic aftereffects were studied. ese pulses are shown in Figure 2 (b) and the measured RMTs are displayed in Figure 3 (c).

Goetz et al. quantified changes in the corticospinal excitability produced by three near-rectangular-shaped cTMS stimuli, one bidirectional and two unidirectional in opposite directions (named RB : Rectangular bidirectional. RU-N : Rectangular unidirectional with initially AP. RU-R : Rectangular unidirectional with initially PA) [28] . ese pulses are shown in Figure 2 (c) and the measured RMTs are displayed in Figure 3 (d).

In [29] , the plasticity effect of unidirectional cTMS pulses in the paired-pulse protocol was investigated. Two types of MT were measured in that research.

#1RMT : resting motor threshold for a positive pulse width of 60 μs and m � 0.2. #2RMT : e pulse magnitude was then set to 120% of #1RMT, which the authors called the "test intensity." ey then measured a threshold pulse width ( PW) that elicits resting MEPs with a 50% probability. According to the results of this study, the average #1RMT is 32% of maximum stimulator output (MSO) and the PW of #2RMT is 47.3 μs, as shown in Figure 4 .

Although we call the pulses by their higher pulse widths in terms of absolute amplitude in the cTMS device, such as 40 μs or 80 μs pulses, in all calculations, the compensatory phases are also included in the waveform and the MT estimation. Selecting the width of the major phase is arbitrary and the opposite phase is selected to reach the almost zero current at the end of each pulse. M ratio has a direct effect on the choice of positive and negative pulse widths to reach this goal. Also, due to some limitations of the cTMS equipment in charging and discharging high voltage DC capacitors, for repetitive or paired pulses protocols, balanced pulses must be defined. e shape of the cTMS pulses was selected based on the experiments defined in our database, and their purpose was to investigate the effect of different pulse shapes on the neural behavior, both monophasic and biphasic pulses. Starting or ending with a compensatory phase was done by researchers in our selected databases, and we have no role in defining pulses. We simulated their selected waveforms and estimated the MT.

e estimated MT values and the experimental values in the defined database can be seen in Figures 3 and 4 . e maximum absolute error between the average of the RMT and the AMT for monophasic stimuli (assumed as the true motor threshold) [5, 27] and the MTestimates made with the HH model was 5% at 30 μs and 4% at 120 μs PW, respectively. e maximum absolute error between the average of the RMT for the three unidirectional stimuli [8] and the estimated MTs was 8% at 40 μs PW. e error between the x + E 2 y + E 2 z ). 4 Computational Intelligence and Neuroscience average of the RMT for one bidirectional and two unidirectional (in opposite directions) stimuli [28] was 3% for the RU-R pulse.

According to the #1RMT protocol in [29] , the 60 μs stimulus intensity is estimated to be 39.2% of MSO to generate an action potential in the HH model. According to the #2RMT protocol, the PW is estimated to be 50 μs to evoke a response. erefore, the maximum absolute error between the #1RMT and #2RMT protocols and the MT estimates obtained with the HH model was 7.2% and 5.7%, respectively. e measurement and estimation results are shown in Figure 4 

In this research, the novel method to predict the RMT and the AMT was introduced. e results of human experiments presented in five previously published studies, which used the cTMS device for their trials, were collected as a database. Based on the outcomes of the computational modeling and [29] with initial AP direction (m � 0.2). e ratio of the peak of the negative phase to the positive phase is called the m ratio. Generally, the stimulus intensity is proportional to the E-field, (a derivative of the B-field), i.e., the current density in the brain and input to the HH model.

Computational Intelligence and Neuroscience the comparison with the database results, it can be concluded that the proposed method can estimate the values of the MT for all three types of monophasic, unidirectional, and bidirectional pulses with different directions and widths. In a similar attempt, Julkunen applied the cumulative distribution function (CDF) with Monte Carlo simulations to determine an MT estimate [9] . Although the precision of that method is higher than in this study (error ≤ 0.9%), its algorithm can only be used for conventional pulses [30] , and it is not possible to predict the MT value by manipulating the stimulus character (such as the waveform, frequency, and directionality).

Recent computational and experimental studies highlight the complexity in the MT estimation with TMS. Wu et al. implemented the cable theory which was adjusted in the detailed compartmental neural cell models [31] . Aberrra et al. utilized the FEM with morphologically realistic models of nerve cells to estimate cortical response to the different magnetic stimuli and quantified the activation threshold for layers 1-6 [32] . Wilson et al. used a neural field approach with a nonlinear model of motoneuron firing to model motor-evoked potentials as a function of amplitude, including identification of motor thresholds [33] . However, despite comprehensive approaches, the increasing sophistication of these computational models dramatically reduces their computational efficiency [34] . In contrast, this present work demonstrates that a straightforward calculation based on the wellestablished HH model and finite-element modeling of electric fields in a simple geometry is sufficient to reproduce experimental results for motor thresholds. However, more complicated approaches may well be required for modeling effects above the threshold, for example, the size of motorevoked potentials. Experimental results in [5, 27] HH simulation results Experimental results in [28] HH simulation results

RU-N RU-R Different pulse shapes (d) Figure 3 : e effect of different pulse shapes on the MT (mean ± standard deviation). Strength-duration curves for (a) RMT measurement results and (b) AMT measurement results for monophasic pulses in [5, 27] (n � 10), which are compared to the simulation results of the HH model, (c) RMT measurement results for unidirectional pulses in [8] (n � 15). (d) Resting motor thresholds for RB, RU-N, and RU-R pulses in [28] (n � 13), in comparison with the modeling results. e stimulus intensity reflects the MT as a percentage of the maximum stimulator output. e maximum voltage of the cTMS device is 2800 V which is equivalent to 100% stimulus intensity. e filled red area shows the standard deviation of the mean, according to the experimental results. All estimated MT values are inside or very close to the area of the standard deviation. 6 Computational Intelligence and Neuroscience

Finally, we emphasize that the model has reproduced experimental results without tuning parameters. We have taken a standard implementation of the HH model and made a physical estimate of the induced current density J stim for a given applied voltage to the coil with no fitting of parameters.

Evidently, the effects of magnetic pulses on neural tissues are complex and are a function of multiple variables, including a contribution from the brain (gyral) geometry [35] , the cortex position [36] , the size of the somatic compartment, and the input resistance [37] . Modeling these parameters requires the implementation of complex equations and possibly FEM for geometric modeling of the gyral crown.

us, the proposed method cannot predict RMT or AMT for individuals but can predict average thresholds for a population of people. It may be possible to individualize the approach by changing the mapping between peak coil voltage and induced current density J stim for different people based on a standardized pulse shape response.

In this study, we focused on the activation kinetics of the voltage-gated ion channels. We have assumed that the input current density to the HH neuron is directly related to the induced electric field due to the TMS pulse. Despite the many simplifications, all of our modeling results have clearly correlated with previous works. Using the assumption that ion channel modeling can predict the overall strength-duration trend with an acceptable error range, it may be possible to estimate the intensity for the stimulus by embedding the differential equations in the hardware.

Additionally, several ions are assumed to be involved in the activation kinetics of the ion channels stimulated by near-rectangular pulses, but a possible connection between the pharmacological impact on motoneuron firing has not been determined yet. Considering that the MEP creation was affected by the administration of neuroactive drugs which are related to neurotransmitters [38] , further pharmacological research is needed to describe the effect of magnetic pulse shapes and the MEP.

For the first time, the suggested approach has been shown to estimate the MT for cTMS stimuli. e MT assessment allows investigators to refine the strengthduration relationship for controllable stimuli. Finding the optimal pulse with the minimum energy required to stimulate the neuron can steer to the plan of a novel magnetic pulse generator equipment with minimal hardware and size.

is issue can be useful in home-TMS equipment, especially now that there is an urgent need for more accessible at-home neurotherapeutics throughout this COVID-19 social isolation period [39] . e prediction of the MT can also minimize the number of experiments required to identify the MT for new magnetic stimuli.

Due to the regulation of the excitability of nerve cells by several ion channels and also membrane capacitance, a critical inquiry for the model accuracy is to find the sensitivity of neural behavior to a variation in HH parameters: which value is most important for the MT estimation? To simplify sensitivity analysis, we focus on RU-N pulse and 7 parameters, which are assumed to change independently and all remaining parameters keep unchanged and ignore changes in equilibrium potentials. Analyzed parameters include equivalence membrane capacitance (C m ), sodium channel conductance (g Na ), potassium channel conductance (g K ), slow potassium channel conductance (g M ), leak channel conductance (g leak ), time constant for the slow potassium channel (τ max ), and TMSinduced current density (J stim ). e calculated results for changing the estimated MT value for this sensitivity analysis are shown in Figure 5 . e range of changes in the parameters has been selected according to the different values measured in different references, biological variation in the parameters [15] , and the possibility of observing the change in the estimated value of the MT. e electrical conductivity variation is assumed to be ±20%, based on the in vitro measurements presented in [40] .

According to the computational sensitivity results, it can be seen that the highest sensitivities are related to the values of membrane capacitance and sodium channel parameters. e sensitivity of the MT to the values of potassium and TMS-induced current density is less than the previous parameters. Also, the HH equations are almost insensitive to fluctuations in leak channel conductance, potassium channel conductance, and time constant for the potassium channel, which are not shown in Figure 5 . Mathematical analysis of the neuronal model, instead of computational analysis, can calculate more accurate results of the sensitivity of the proposed method [15, 41] .

cTMS equipment has opened up novel opportunities to investigate the effect of magnetic pulse shape on the human nervous system. In this study, we introduce a novel method for determining MTs for cTMS device pulses based on the HH model, as mathematical descriptions of neuronal behavior that is exposed to the magnetic stimulation. e HH model was selected as an interconnection of a biochemical network with mathematical neuroscience. e results of this computational study highlighted a strong contribution to the activation kinetics of the voltage-gated ion channels for the motor threshold induced by TMS. Validation with previously published experimental threshold data confirms that the proposed method captures the measured population-averaged MT value at high precision, without the need to perform finite-element modeling in realistic head geometries or consideration of large networks or populations of cells. e proposed prediction method will allow a rapid estimation of the motor threshold to aid decisionmaking by researchers and clinicians to determine the motor threshold and improve the safety implications of preintervention probing. Furthermore, the proposed mathematical framework might be a promising tool to explore efficient magnetic pulses to design a home-TMS machine.

No data were used in the study.

e authors declared that they have no conflicts of interest to this work.

Non-invasive electrical and magnetic stimulation of the brain, spinal cord, roots and peripheral nerves: basic principles and procedures for routine clinical and research application. An updated report from an I.F.C.N. Committee

TMS brain mapping of the pharyngeal cortical representation in healthy subjects

Deep-brain transcranial stimulation: a novel approach for high 3-D resolution

Numerical modeling of plasticity induced by quadri-pulse stimulation

Controllable pulse parameter transcranial magnetic stimulator with enhanced circuit topology and pulse shaping

Programmable transcranial magnetic stimulation-a modulation approach for the generation of controllable magnetic stimuli

Paralleling insulated-gate bipolar transistors in the H-bridge structure to reduce current stress

Neuronal tuning: selective targeting of neuronal populations via manipulation of pulse width and directionality

Mobile application for adaptive threshold hunting in transcranial magnetic stimulation

Active and resting motor threshold are efficiently obtained with adaptive threshold hunting

Controllable pulse parameter transcranial magnetic stimulator with enhanced circuit topology and pulse shaping

eoretical analysis of transcranial magneto-acoustical stimulation with hodgkin-huxley neuron model

A modified Hodgkin-Huxley model to show the effect of motor cortex stimulation on the trigeminal neuralgia network

Dynamic redistribution of plasticity in a cerebellar spiking neural network reproducing an associative learning task perturbed by TMS

Cellular function given parametric variation in the Hodgkin and Huxley model of excitability

Parameter estimation for bursting neural models

Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons

e New York Head-A precise standardized volume conductor model for EEG source localization and tES targeting

Electric field depth-focality tradeoff in transcranial magnetic stimulation: simulation comparison of 50 coil designs

Measurement of transcranial magnetic stimulation resolution in 3-D spaces

Fast estimation of transcranial magnetic stimulation motor threshold

TMS motor mapping: comparing the absolute reliability of digital reconstruction methods to the golden standard

Age and sex differences in human motor cortex input-output characteristics

Chapter 23 pharmacology of TMS

e cortical motor threshold reflects microstructural properties of cerebral white matter

Half sine, monophasic and biphasic transcranial magnetic stimulation of the human motor cortex

Effect of coil orientation on strength-duration time constant and I-wave activation with controllable pulse parameter transcranial magnetic stimulation

Enhancement of neuromodulation with novel pulse shapes generated by controllable pulse parameter transcranial magnetic stimulation

Strength-duration relationship in paired-pulse transcranial magnetic stimulation (TMS) and its implications for repetitive TMS

Chapter 2 TMS and threshold hunting

Cortical neuron activation induced by electromagnetic stimulation: a quantitative analysis via modelling and simulation

Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons

Modeling motor-evoked potentials from neural field simulations of transcranial magnetic stimulation

Which model to use for cortical spiking neurons

Impact of the gyral geometry on the electric field induced by transcranial magnetic stimulationfield induced by transcranial

Where does TMS stimulate the motor cortex? combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position

Patch-clamp recordings of rat neurons from acute brain slices of the somatosensory cortex during magnetic stimulation

e pharmacology of neuroplasticity induced by non-invasive brain stimulation: building models for the clinical use of CNS active drugs

Treating the mental health effects of COVID-19: the need for at-home neurotherapeutics is now

In vitro electrical conductivity of seizing and non-seizing mouse brain slices at 10 kHz

Neuronal behaviors: a control perspective

Acknowledgments is work was supported by a program grant from the MRC Brain Network Dynamics Unit at Oxford (MRC MC_UU_0003/3) and supplemental funding to Denison by the Royal Academy of Engineering and Wendt by an MRC iCASE Fellowship. Figure 1S : diagram of key circuit elements of the controllable-TMS device.