Thermal effects on neurons during stimulation of the brain

We observed the activities of cortical layer V pyramidal cells from male Sprague-Dawley rat brain slices in vitro, using a loose-patch attachment [37] to reduce the possibility of electrical shunting into the neuron through the electrode during magnetic stimulation as in [33]. Stimulation was applied by driving micromagnetic coils with either AC or DC current adjusted for equivalent power dissipation ($P = I_\mathrm{RMS}^2R$) and thus equivalent temperature change measured near the patched cell (figures 1(a) and (b) see Methods). Unexpectedly, we observed that the pattern of spike rate changes with AC or DC stimulation could be identical, with transient suppression at the onset of stimulation current and transient hyperexcitability after removal of current, accompanied by similar changes in temperature at the level of the chamber floor near the coil and cells (figures 1(c) and (d)).

Zoom In Zoom Out Reset image size

We reduced the amplitude of the stimulation current in the micromagnetic coils to limit the temperature increase near the patched cell to about 1 ^∘C for 180 s, and then allowed adiabatic cooling back to the 30 ^∘C bath temperature baseline upon cessation of coil stimulation (figure 2(d)). We verified this equivalent temperature increase by placing a temperature sensor near the patched cell (figure 1(a)), and found the time constant of temperature increase and decrease to be on the order of 5 s (figure 2(d)).

Zoom In Zoom Out Reset image size

During AC stimulation the coil was driven with continuous sinusoids to induce changing magnetic fields at 50 Hz, 500 Hz, 5 kHz and 200 kHz, while controlling the current amplitude to maintain the same power dissipation and 1 ^∘C temperature increase across all tested frequencies. In order to isolate thermal effects from magnetic stimulation, we delivered AC and DC current delivered in a block design, alternating AC-DC-AC with DC-AC-DC stimulation blocks with flanking control blocks without stimulation (figure 2(b)). Each stimulation block was subdivided into 10 s sub-blocks (per 60 s block). This enabled us to determine whether any magnetically driven electrical current induction affected the cell's action potential frequency or height within each sub-block.

We identified a statistically significant transient suppression of spikes immediately following the onset of stimulation (mean spike density $0.484\pm0.294$ normalized to baseline), and transient hyperactivity following cessation of stimulation (mean spike density $1.295\pm0.860$ normalized to baseline, figures 2(a) and (c)). This effect was independent of whether AC or DC stimulation was applied to the coils, and the equivalent power dissipation was reflected in identical temperature increase and decrease profiles (figure 2(d), see also figure A1 in the appendix).

Across all AC stimulation frequencies, we found no statistically significant differences in spike height or spike frequency within a given sub-block (table A2). That is, the changes seen at the onset and offset of stimulation appeared independent of driving frequency applied.

We next determined whether spike timing was changed by the AC stimulation sinusoids. This was done by testing whether the spikes tended to occur at particular phases of the AC sinusoids. The statistical distributions of the phases of the spike times were evaluated using the Raleigh Z statistic for circular distributions (see Methods). No AC blocks showed statistically significant spike phase preference compared to the DC blocks (figure A7). We also measured the directionality of spike phases in 3 s time windows to test whether spike entrainment might occur more often in AC blocks than in the DC blocks over short periods of time. To reduce the chance of type II error, we identified candidate time windows showing entrainment by comparing data with randomly shuffled surrogate interspike intervals drawn from the same data blocks and analyzing the spectral power density of the ensemble (see Methods). The populations of Rayleigh Z p-values calculated from these candidate time windows for AC and DC blocks were compared to each other using the Anderson Darling test [38] and Wilcoxon ranksum test. The results (table 1) showed that there was no statistically significant difference between the spikes in DC and all frequencies of AC blocks. This extensive set of the statistical analysis of experiments and controls was required to prove that the effects on neural activity were due to thermal effects rather than any detectable trace of magnetic induction of electrical currents interacting with the neurons.

At a distance of 150 µm from the coil, the magnitude of the induced electric field strength in the tissue ranged from 0.08 µV mm⁻¹ to 0.3 mV mm⁻¹ for 50 Hz–200 kHz respectively (see Methods and figure A2(b). The experimental [39] and theoretical [40] limits for field interactions with the spike generating mechanism of neurons is about 0.1 mV mm⁻¹ at driving rates of stimulation below 50 Hz. Entrainment of single neurons and networks is readily observed at these lower frequencies [39, 41]. In the present analysis, the lack of detection of any interaction of the neuronal spiking with driving fields is consistent with the subthreshold fields generated at the lower frequencies tested.

In addition to the thermal transients due to conduction of coil Joule heating to the tissue, the SAR in the tissue under the experimental conditions presented here requires consideration. The SAR is computed as

$$\begin{align} \mathrm{SAR} = \frac{\sigma(r)|E(r)|^2}{2\rho(r)}, \end{align} \tag{ 1 }$$

where σ is the tissue electrical conductivity, E is the induced electric field in the tissue and ρ is the tissue density [42]. We used tissue conductivity and density values reported in [42] as 0.276 S m⁻¹ and 1035.5 kg m⁻³ respectively. The induced electric field was estimated from a COMSOL finite element frequency quasi-static analysis for the range of input frequencies tested experimentally (see Methods, figure A2 and table 3). The computed range of the maximum local SAR in the tissue from the COMSOL model was $4.50\times10^{-11}$–$7.64\times10^{-4}$ W kg⁻¹, which was negligible compared with the established allowable localized SAR safety limit in the head of 2–4 W kg⁻¹ [42]. This computation, taken in combination with our results above, confirms that in the thermal profile of the stimulation Joule heating dominates the effects on the tissue over the SAR from the induced field absorption.

To understand the thermal effects observed in our experiments we modified a computational model based on [43], which extends the original Hodgkin-Huxley formalism [44] to include relevant structural micro-anatomy, conservation of charge and mass, energy balance, and volume changes. The model consists of a single neuronal compartment surrounded by a thin extracellular space that is connected to the bath in vitro (or to a capillary in vivo) through ionic diffusion (figure A3). These biophysical components are required to unify a range of phenomena such as spikes, seizures and spreading depression, which the original Hodgkin-Huxley framework does not characterize well for mammalian neurons. The Nernst potentials that underlie the transmembrane voltage dependency on ionic concentration differences are proportional to temperature. In addition to the thermal effect inherent in the Nernst potentials, we identified eleven membrane processes in the model that could be important in explaining the thermal effects observed: maximal Na-K pump rate, maximal Na⁺ and K⁺ channel conductances, rate of gate dynamics ($\mathrm{d}n/\mathrm{d}t$, $\mathrm{d}m/\mathrm{d}t$, $\mathrm{d}h/\mathrm{d}t$), membrane leak conductances for Na⁺, K⁺ and Cl⁻, and strength of cotransporters NKCC1 and KCC2. The effects of these thermally sensitive processes on membrane potential can be categorized as hyperpolarizing (pump strength, K⁺ and Cl⁻ conductances, and Nernst potentials), depolarizing (Na⁺ conductances), or charge-neutral (gating kinetics and cotransporter strengths).

Although biological processes have a complex relationship to temperature [45], individual cellular processes generally accelerate when temperature is raised a small amount. We used Q₁₀ factors to model these accelerated rates. The Q₁₀ factor of a process y, describes the change of reaction rate when the temperature is raised by 10 ^∘C. For an arbitrary temperature change, $T-T_0$, the new maximal rate of a process y(T), with initial process rate $y(T_0)$ is calculated as:

$$\begin{align} y(T) = y(T_0)Q_{10y}^{(T-T_{0})/(10\ ^{\circ}\mathrm{C})}. \end{align} \tag{ 2 }$$

When all temperature dependencies in the model are used to simulate the effect of temperature increase or decrease from baseline, transient spike suppression or hyperactivity is observed as in the experiments (figures 3(a)–(d)). We next examined the effect of individual model components on cellular activity. Increasing the temperature-mediated effects of the hyperpolarizing processes transiently silences the cell, while increasing the effects of the depolarizing processes induces hyperactivity (figure A4). Decreasing the temperature-mediated effects of these parameters has opposite results, while changing charge-neutral parameters does not create significant transient behaviors. Since our in vitro data show silencing at the application of heat and hyperactivity at the removal (figure 2(a)), the hyperpolarizing components contained in the model appear to dominate the cell's thermal response.

Zoom In Zoom Out Reset image size

We next use the simplifying assumption that similar membrane process categories (ion transport, Nernst potentials, etc) will have similar Q₁₀ factors (table 5). When the Q₁₀ factors are the same, the cellular process rates will change by the same proportion for any change in temperature. We simulated a temperature rise as a step increase of temperature affecting each grouped category of membrane processes (hyperpolarizing or depolarizing) by 105%, 110%, and 140% followed by a return to baseline values (figure A5). Changing either the leak conductances or the voltage-dependent channel conductances induced hyperactivity at the onset of heating, rather than silencing as observed in the experiments. Thus, under assumptions of similar Q₁₀ factors, the depolarizing Na⁺ currents would dominate, and this is inconsistent with the experimental findings. Relaxing the Q₁₀ similarity constraint for similar biophysical processes, such as leak currents, still reveals that $G_\mathrm{NaLeak}$ would dominate the hyperpolarizing effects of $G_\mathrm{KLeak}$ and $G_\mathrm{ClLeak}$ for similar temperature increases (figure A6).

The voltage gated channel conductances, $G_\mathrm{K}$ and $G_\mathrm{Na}$, do not substantially suppress the cell activity during heating even when their values are set to maximal (figure A4). In contrast, the Na-K pump current demonstrates transient suppression with a temperature increase, and transient hyperactivity with temperature return to baseline (figure A4). In addition, we found that increasing the temperature in the Nernst potential equations for Na, K, and Cl gave a similar transient suppression and hyperactivity as the Na-K pump (figure A4). The charge-neutral co-transporters had no effect on model cellular activity.

When the model is periodically firing, without thermal perturbation, it is in a stable limit cycle. This model contains variables characterized by fast and slow time scales (figures 3(a)–(l)). The fast variables are the membrane voltage V and the three gating variables m, n, and h. The slow variables are the analyte concentrations inside and outside the cell. These concentrations define the environment in which action potentials are generated, establish the Nernst potentials and set the pump currents, thereby modifying the spiking behavior. When a dynamical system in a stable limit cycle is perturbed by an abrupt change of parameters, such as a step change in temperature, the system must then traverse through the state space from the old stable state to the new one (provided the new parameters also define a new stable basin of attraction).

There are substantial qualitative differences in spiking activity when temperature is changed slow or fast. The spiking behavior of the model in figures 4(a) and (b) demonstrates the effect of the speed of the temperature change on the model cell's spiking. Faster temperature time constants create transient silencing and hyperactivity in the model analogous to the temperature effect on in vitro neurons in figures 1(c), (d) and 2(a). In figures 4(c) and (d), the model trajectories resulting from faster rates of temperature changes take large detours in state space, rather than the more direct paths of slower changes. Slow changes in temperature do not alter the equilibrium of the system much, and allow the faster variables to come to equilibrium at all times adiabatically. Slow changes in temperature can even preserve the stable limit cycle, shown as broad ribbons within the trajectories in state space reflecting the oscillations of the spiking limit cycles (time constants of 50 and 500 s in figure 4). With slower temperature changes the model does not display transient silencing or hyperactivity. But for fast temperature changes, the new temperature driven equilibrium can be far away from the present state, and needs to be reached as the system readjusts by its own internal fast and slow dynamics. Such substantial trajectory excursions through state space, arriving at the same final equilibrium point, may take the system far from its original limit cycle frequency and may abolish the spiking entirely as shown in the model and experiment.

Zoom In Zoom Out Reset image size

These transient effects can be described by a type of bifurcation whose stability is characterized as a saddle-node. First, we separate the model state vector into fast and slow variables. Then we treat the trajectory of the slow variable state vector $\pmb{C}(t_\mathrm{s}) = [[\mathrm{Na}^+]_\mathrm{i},[\mathrm{Na}^+]_\mathrm{o},[\mathrm{K}^+]_\mathrm{i},[\mathrm{K}^+]_\mathrm{o},[\mathrm{Cl}^-]_\mathrm{i},[\mathrm{Cl}^-]_\mathrm{o},[\mathrm{O}_2]]$ as a function of slow time $t_\mathrm{s}$. We also assume that the gating variables with their fast time scales are at their steady state values. Then, we plot $\dot{V}(V,\pmb{C}(t_\mathrm{s}))$ vs V (figures 3(i)–(l)). The x-intercepts of this graph are the equilibria of V. We can understand how temperature affects the spiking behavior of the model by studying these fixed points. At the onset of a modeled temperature increase (figures 3(i) and (j)), there exists one stable equilibrium point (P1) at −63.0 mV and two unstable equilibria points (P2 and P3) at −58.9 mV and −41.1 mV (see figure A8). The stable fixed point at P1 corresponds to quiescent cell behavior as seen in the in vitro patch-clamp data in figures 2(a) and (b), where the cell ceases firing upon warming.

We determine the stability of the gating variables at the equilibria by calculating the eigenvalues of the Jacobian matrix of V and the gating variables m, h, and n at each equilibrium. The Jacobian and its components are shown in appendix 'Stability of equilibria in the $\dot{V}$ curve'. All eigenvalues must have negative real parts for the equilibirum to be stable. This analysis shows that P1 is the only stable equlibirum.

We can verify this result near these equilibrium voltages, where we know that $\mathrm{d}m/\mathrm{d}t \gg \mathrm{d}n/\mathrm{d}t \gg \mathrm{d}h/\mathrm{d}t$ (figure A9). To illustrate the stability of the equilibrium points, we generated field maps for each of the gating variables versus V, assuming that for each gating variable, faster gating variables are at their steady state values for the corresponding voltage ($x_{\infty} (V) = \alpha_{x}(1-x)-\beta_{x}$, where x = m, n, h), and slower gating variables are at their steady state values for the equilibrium voltage ($x_{\infty} (V_\mathrm{eq})$) in figure A10 [46]. For example, in calculating the field map of n vs V, we assumed that m will be at $m_{\infty} (V\,)$, and h will be at $h_{\infty} (V_\mathrm{eq})$. This analysis confirmed that P1 is the only stable equilibirum.

Regarding the slow variables, we treat the slow time, $t_\mathrm{s}$, as the bifurcation parameter for the equilibria. As $t_\mathrm{s}$ progresses following warming, and by extension the analyte concentrations shift, the $\dot{V}$ curve moves up. Eventually the two equilibria (points P1 and P2) collide, creating a saddle-node bifurcation. Once the stable equilibrium disappears (the quiescent state), the cell is able to depolarize further to initiate action potentials. In a similiar manner the $\dot{V}$ curve can also be used to predict hyperactivity upon cooling. Figures 3(k) and (l) show the $\dot{V}$ curve when the heat is removed. The distance between the local minimum of the curve and the line $\frac{d}{\mathrm{d}V}(\dot{V\,}) = 0$ is representative of the cell's depolarization speed, and consequently, the frequency of spikes. When the curve is below this line, the distance between the curve minimum and $\frac{d}{\mathrm{d}V}(\dot{V}) = 0$ represents how resistant the cell is to spiking.

To further characterize the biophysics, we reduced the complexity of the ion species relationships (see Methods) and used numerical continuation analysis to explore bifurcations at a broader range of temperatures. While models of small changes in temperature showed transient changes but not qualitative differences in longer term equilibrium spiking behavior, bifurcation analysis of the reduced model at higher temperature revealed a subcritical Hopf bifurcation around 37.26 ^∘C (figure 5(a)). Creation of this bifurcation only requires the Nernst potential to change with temperature; it exists even when all other Q₁₀'s are set to 1. When the temperature dependence of the Nernst potential is introduced to the basic HH model without any charge or ion conservation, the Hopf bifurcation point changes to a saddle-node, but the stable branch at high temperature still emerged.

Zoom In Zoom Out Reset image size

While all membrane parameters affect the location of this bifurcation point, the KCC2 cotransporter expression level also qualitatively changes the behavior of the cell as the temperature is increased. In the less mature developing brain, it is known that KCC2 levels are not as high as in the mature brain, and febrile seizures are more commonly observed with hyperthermia (typically at the onset of fever as temperature is increasing). In the full model with normal mature KCC2 levels (100% strength in our model), spiking is suppressed above 36 ^∘C. At normal KCC2 levels, the frequency and temperature describe an inverse relationship in the full model, while at lower KCC2 levels, the relationship becomes positively correlated and spiking is not suppressed (figure 5(b)). Furthermore, the equilibrium chloride level, $E_\mathrm{Cl}$, becomes elevated (more depolarized) at lower levels of KCC2 expression, reducing the effectiveness of inhibitory transmitter on the cell (figure 5(c)). These effects of KCC2 combine to render spiking at elevated temperatures more prominent in the immature brain with lower levels of KCC2 expression.

We performed a series of experiments delivering thermal and magnetic stimulation to neocortical tissues with custom microscale electromagnetic coils. The experimental protocol was designed from three-epoch trials alternating from DC current (thermal stimulation) to AC current (magnetic plus thermal stimulation) and reversed as shown in figure 2(b) (see also figure A1). This stimulation technique was implemented to control for the thermal energy delivered to the tissue via Joule heating of the micro-coil during the brief but high intensity current pulses required to induce magnetic and electric fields. The AC stimulation amplitude was tuned to produce the same measured temperature as the DC current near the cell, thus eliminating the temperature step that would have been otherwise coincident with the onset of magnetic fields produced by the AC currents. After the block of AC current stimulation, we switched back to DC current to remove the cooling temperature step from coinciding with the offset of AC stimulation. In order to avoid any ordering bias in our experimental design we also swapped the order of AC and DC blocks.

All experiments were performed in accordance with and approval from the Institutional Animal Care and Use Committee of The Pennsylvania State University. Neocortical slices were obtained from male Sprague-Dawley rats aged P11 to P14. Briefly, the animals were deeply anesthetized with diethyl-ether and decapitated, the brain removed and coronal slices were sectioned from occipital neocortex.

After decapitation, the whole brain was quickly and carefully removed to chilled (4^∘C) artificial cerebrospinal fluid (ACSF) for 60 s containing the following (in mM): 126 NaCl, 2.5 KCl, 2.4 CaCl₂, 10 MgCl₂, 1.2 NaH₂PO₄, 18 NaHCO₃, and 10 dextrose. ACSF was saturated with 95% O₂ and 5% CO₂ at room temperature for 1 h before the dissection with an osmolality range from 295 to 310 mOsm and pH from 7.20 to 7.40. Coronal occipital cortical slices were cut with a vibratome on the rostro-caudal and medio-lateral coordinates of bregma −2 to −8 mm and lateral 1–6 mm, respectively. The first cut was made 1100 µm deep from the caudal surface and discarded, and three slices were taken from the brain each 300 µm thick. After cutting, slices were transferred to a chamber containing ACSF and allowed to recover for 30 min at 32 ^∘C–34 ^∘C, then incubated at room temperature (20 ^∘C–22 ^∘C) for an additional 30 min prior to recording. In order to examine the modulation of spiking in layer 5 pyramidal neurons, network spiking activity was supported by modified ACSF consisting of (in mM) 121.25 NaCl, 6.25 KCl, 1.5 CaCl₂, 0.5 MgCl₂, 1.25 NaH₂PO₄, 25 NaHCO₃, and 25 dextrose, perfused over the slice at a rate of 2.8–3.0 ml min⁻¹ at 30 ^∘C during recordings.

The micro-coils (figure 1(b)) were fabricated on a 150 µm thick glass cover slip that formed the floor of our slice recording chamber. First the cover slip is cleaned with acetone, isopropyl alcohol and rinsed with deionized water followed by an O₂ plasma ashing to remove organic residue. Next, a conductive seed layer of 20 nm thick Cr followed by 120 nm Au is evaporated for copper electroplating. To enable electroplating of thick copper coils a 28 µm thick MEGAPOSIT SPR 220 photoresist (Shipley Company, Marlborough, MA) is spun on the slides in two coats at 1500 RPM followed by a softbake. To minimize cracking and bubbling of the resist the softbake is conducted on two hotplates first at 95 ^∘C for 90 s, then 115 ^∘C for 120 s and finally back down to 95 ^∘C for 90 s. The photoresist is exposed at a dose of 1200 mJ cm⁻² on a SÜSS MA/BA6 contact aligner (SÜSS MicroTec SE, Garching, Germany) in 15 cycles with 15 s delays to reduce the risk of bubbling and cracking of the resist. This is followed by a minimum 30 min wait time before being developed in MICROPOSIT MF CD-26 (Shipley Company, Marlborough, MA) for 11 min. A final O₂ plasma ashing is then performed to remove any remaining photoresist in the coil traces. The copper coil traces are then electroplated in a beaker using a high purity 99.99% copper anode and a Technic Copper FB Bath plating solution (Technic Inc. Cranston, RI) at 7 mA with a 20% duty cycle for 4 h, followed by palladium electroplating to cap the copper with a non-oxidizing noble metal. The palladium is electroplated using Technic Pallaspeed RTU bathing solution (Technic Inc. Cranston, RI) using a Pt anode (Pd is dissolved in the solution) at 7 mA with a 20% duty cycle for 20 min. The photoresist is stripped using RemoverPG (MicroChem, Newton, MA) and the thin Cr/Au seed layer is etched away using an Ar/SF6 reactive ion etch. Finally an electroless Au is deposited using a gold potassium cyanide solution at 90 ^∘C for 2 h. The coils are measured to be between 23 µm and 26 µm thick. This type of non-uniformity is typical of a beaker electroplating set-up. The resistance of the micro-coils used were measured between 0.43 and 0.60 Ω with a multimeter (Keithley 2000, Keithley Instruments, Solon, OH) before and after each experimental trial.

The recording chamber (RC-22C; Warner Instruments, Holliston, MA) was fitted with a 150 µm thick cover slip and the chamber floor modified with lithographically deposited micro coils on the bottom (dry) side of the glass as shown in figures 1(a) and (b). Patch pipettes (4–6 MΩ) were pulled from thick-walled borosilicate glass capillaries with filaments (OD = 1.5 mm, ID = 0.86 mm; Sutter Instruments, Novato, CA) and filled with perfusate. Layer V pyramidal cells were targeted under visual control (figure 1(a)) and the patch pipette was adhered to the soma of the cell in loose patch configuration [37] (5–10 MΩ seal) in order to minimize the possibility of shunt current entering the cell from the pipette during magnetic stimulation. To drive magnetic stimulation, a function generator (Keithley 3390; Keithley Instruments, Solon, OH) supplied the AC input signal to an audio amplifier (PB717X, gain 2.16X; Pyramid Inc. or AE Techron 7224, gain 1.43X; AE Techron, Elkhart, IN, for high frequency stimulation up to 200 kHz) which then amplified the signal sent to the coil. For thermal-only stimulation a DC power supply (Keithley 4590; Keithley Instruments, Solon, OH) was connected to the micro coil to supply heat to the tissue (stimulation details found in table 2). Neuronal spiking was recorded (MultiClamp 700A; Axon Instruments, San Jose, CA) under perfusion with the modified ACSF and only actively spiking cells were tested with experimental stimulation. A Labview interface was used to low-pass filter and amplify the thermal signature recorded: (a) at the surface of the chamber floor glass (wetted side) to measure the maximum temperature experienced by the tissue, and (b) near (∼50 µm) the patched cell using multiple highly sensitive (68 µV ^∘C⁻¹, diameter 76 µm) micro E-type thermocouples (5SRTC-TT-E-40-36, Omega Engineering Inc.) as shown in figure 1(a).

Table 2. Summary of stimulation currents used during each frequency tested. $N_\mathrm{animals}$ and $N_\mathrm{cells}$ are the number of animals and patched cells for each frequency tested. The micro-coils used during stimulation ranged from 0.43 to 0.60 Ω.

Experimental stimulation parameters
Frequency	$N_\mathrm{animals}$	$N_\mathrm{cells}$	Coil Current (DC)	Input Voltage (AC)
50 Hz	5	12	0.575 ± 0.019 A	0.417 ± 0.012 V
500 Hz	8	20	0.571 ± 0.015 A	0.415 ± 0.016 V
5 kHz	7	24	0.565 ± 0.020 A	0.417 ± 0.025 V
200 kHz	9	16	0.751 ± 0.093 A	0.701 ± 0.077 V

Figure A1 shows a typical experimental result. During the blue colored control blocks A and E, no current was driven through the micro-coils, resulting in no additional heat or electromagnetic fields. In the green and red blocks B though D, the micro-coil is being driven in either AC or DC, producing heat and electromagnetic fields when AC driven. The voltage placed across the micro-coil is shown in the top panel of figure A1(a), and the temperature measured at the probe near the patched cell is shown in the middle panel in figure A1(a). The AC and DC driving stimuli were both voltage adjusted until the temperature measured near the patched cell was ∼1^∘C higher than baseline for each stimulus frequency tested (50–200 kHz).

Spike entrainment refers to the phenomenon where spikes preferentially fire at certain phases of the applied AC sinusoidal stimulation. We used the Rayleigh Z test as the primary statistical tool to detect spike entrainment. For a set of n phase measurements, the Rayliegh Z metric is calculated by:

$$\begin{align} z = n\mathbf{r}^2 , \end{align} \tag{ 3 }$$

where r is a 2 dimensional vector with X and Y components,

$$\begin{align} \mathbf{X}& = \frac{\mathop\sum\limits_{i = 1}^{n} \cos(a_i)}{n} \end{align} \tag{ 4 }$$

$$\begin{align} \mathbf{Y} & = \frac{\mathop\sum\limits_{i = 1}^{n} \sin(a_i)}{n} . \end{align} \tag{ 5 }$$

The null hypothesis of Rayleigh Z testing is that there is no preferred phase. The p-value can be approximated as $\exp( \sqrt{1+4*n+4*(n^2-R_z *n)}$ $-(1+2*n) )$, where R_z is the Rayleigh Z statistic, and n is the number of measurements [84]. We can also calculate Rayliegh Z for the DC blocks, using pseudophases. Pseudophases are measured by overlaying an AC block sinusoid on the DC spike data, and measuring the spike phases. These pseudophases serve as a control for the phases measured in the AC blocks.

First, we tested whether any AC blocks showed a preferential phase. Using the Benjamini–Hochberg procedure to control the false discovery rate at 10%, we found no blocks had statistically significant preferred phases.

We also accounted for the possibility that spike phase entrainment exists for short periods of time. For each block of phases and pseudophases, 99 surrogate spike data were generated by randomly shuffling the inter-spike intervals. For the ensemble of real and surrogate data, the time dependent power densities at the stimulation frequency in a moving window of 3 s were calculated using the Chronux package [85] (http://chronux.org). The confidence limit was set at 95th percentile of the ensemble, and we noted the temporal locations in each AC block where the power density of the real data passed the confidence limit. We grouped the time windows at these temporal locations with others if any part of them overlapped. We then obtained groups of p-values from Rayleigh-Z tests performed on the resulting combined candidate time windows, and again controlled for false discovery rate by applying the Benjamini–Hochberg method separately for phases and pseudophases at each stimulation frequency. If this returned any true-phase candidates as potentially significant at a particular stimulation frequency, we compared the distribution of the Rayliegh-Z p-values calculated from phases to those calculated from the pseudophases at that frequency in two ways: First, the Anderson-Darling test was used to test if the two samples come from different distributions. Second, a one-tailed Wilcoxon rank-sum test was used to see if the two distribution have different medians. Both tests failed to show greater degree of spike entrainment in AC blocks compared to DC blocks. Results are summarized in table 1.

We performed a finite element simulation of the induced electric and magnetic fields due to sinusoidal currents driven through our coil in COMSOL multiphysics (COMSOL Inc., Los Angeles, CA). In the simulation the magnetic field (figure A2(a)) due to sinusoidal current flow in the coil and the induced electric field (figure A2(b)) in the surrounding media were computed. We then computed the range of SAR values (equation (1)) occurring inside the tissue domain for the stimulation parameters we used in our in vitro experiments. The power dissipated in the coil was evaluated in the model to be 87 mW under our experimental conditions. The model consists of a 2D axis-symmetric revolved geometry comprising the coil trace copper conductor surrounded by air, glass cover slip and the brain tissue (1 mm × 1 mm total axis-symmetric cross-section) as shown in figure A2(a). A frequency-domain quasi-static approximation was used in the computation, along with a mesh refinement to obtain adequate solution accuracy and spatial resolution. All model parameters used in the simulation were taken from tables 2 and 3.

Table 3. Summary of material properties used in the finite element computation of electromagnetic field induction.

Material properties of the brain tissue and coil
Property	Copper	Glass	Air	Gray matter [42]
σ-electric conductivity (S m−1)	5.998 × 107	1 × 10−14	0	0.276
-relative permittivity	1	4.2	1	12 000
ρ-density (kg m−3)	8960	2210	1.204 × 10−3	1035.5

Our model is based on the unification model from Wei et al 2014 [43]: a single compartment spherical neuron model that expands the Hodgkin-Huxley model [44] to account for relevant structural micro-anatomy, conservation of charge and mass, energy balance, and volume changes. It consists of a single compartment spherical neuron with a thin extracellular space that is connected to the bath in vitro (or to a capillary in vivo) through a diffusion channel.

The membrane potential, V, is defined by:

$$\begin{align} C\frac{\mathrm{d}V}{\mathrm{d}t} &= -I_\mathrm{Na}-I_\mathrm{K}-I_\mathrm{Cl}-I_\mathrm{pump}/\gamma \nonumber \\ I_\mathrm{Na} &= G_\mathrm{Na}m^3h(V-E_\mathrm{Na})+G_\mathrm{NaLeak}(V-E_\mathrm{Na}) \nonumber \\ I_\mathrm{K} &= G_\mathrm{K}n^4(V-E_\mathrm{K})+G_\mathrm{KLeak}(V-E_\mathrm{K}) \nonumber \\ I_\mathrm{Cl} &= G_\mathrm{ClLeak}(V-E_\mathrm{Cl})\nonumber \\ I_\mathrm{pump} &= \rho_\mathrm{max}*\frac{1}{1+\mathrm{exp}((25-[\mathrm{Na}^{+}]_\mathrm{i})/3}\nonumber \\ &\quad *\frac{1}{1+\mathrm{exp}(3.5-[\mathrm{K}^{+}]_\mathrm{o})}\nonumber \\ &\quad *\frac{1}{1+\mathrm{exp}((20-[\mathrm{O}_2])/3} \nonumber \\ \frac{\mathrm{d}q}{\mathrm{d}t} &= \alpha_{q}(1-q)-\beta_{q}q,\;\;\;\;q = m,h,n \end{align} \tag{ 6 }$$

where $I_\mathrm{Na}$, $I_\mathrm{K}$, and $I_\mathrm{Cl}$ are the ohmic sodium, potassium, and chloride currents, respectively. $I_\mathrm{Na}$ and $I_\mathrm{K}$ consist of voltage-gated and leak currents, while $I_\mathrm{Cl}$ consists solely of leak current. $G_\mathrm{Na}$ and $G_\mathrm{K}$ are maximal Na⁺ and K⁺ voltage-gated conductances, and $G_\mathrm{NaLeak}$, $G_\mathrm{KLeak}$, and $G_\mathrm{ClLeak}$ are Na⁺, K⁺, and Cl⁻ leak conductances, respectively. We use conversion factor $\gamma = S/(Fv_\mathrm{i})$ to convert $I_\mathrm{pump}$ from (mM s⁻¹) to (µA cm⁻²), where S is the surface area of the cell, F is the Faraday constant, and $v_\mathrm{i}$ is the intracellular volume of the cell.

Reversal potentials are calculated with Nernst equations:

$$\begin{align} &E_\mathrm{Na} = \frac{\mathrm{RT}}{zF}\ln \left( \frac{[\mathrm{Na}^{+}]_\mathrm{o}}{[\mathrm{Na}^{+}]_\mathrm{i}} \right)\nonumber\\ &E_\mathrm{K} = \frac{\mathrm{RT}}{zF}\ln \left( \frac{[\mathrm{K}^{+}]_\mathrm{o}}{[\mathrm{K}^{+}]_\mathrm{i}} \right)\nonumber\\ &E_\mathrm{Cl} = \frac{\mathrm{RT}}{zF}\ln \left( \frac{[\mathrm{Cl}^{+}]_\mathrm{i}}{[\mathrm{Cl}^{+}]_\mathrm{o}} \right). \end{align} \tag{ 7 }$$

The KCC2 cotransporter was modeled as:

$$\begin{align} I_\mathrm{KCC2} = U_\mathrm{KCC2} \ln \left( \frac{[\mathrm{K}^+]_\mathrm{i}[\mathrm{Cl}^-]_\mathrm{i}}{[\mathrm{K}^+]_\mathrm{o}[\mathrm{Cl}^-]_\mathrm{o}} \right). \end{align} \tag{ 8 }$$

Values and descriptions of the parameters are given in table 4, and the full description of the model is given in the appendix.

Table 4. Values of parameters used in the model.

Parameter	Base value	Description
T	30 ∘C	Temperature
C	1 µF cm−2	Membrane capacitance
$G_\mathrm{Na}$	45 mS cm−2	Maximal conductance of sodium current
$G_\mathrm{K}$	25 mS cm−1	Maximal conductance of potassium current
$G_\mathrm{NaLeak}$	0.0247 mS cm−2	Conductance of leak sodium current
$G_\mathrm{KLeak}$	0.05 mS cm−2	Conductance of leak potassium current
$G_\mathrm{ClLeak}$	0.1 mS cm−2	Conductance of leak chloride current
β0	7	Ratio of the initial intra-/extracellular volume
$\rho_\mathrm{max}$	0.8 mM s−1	Maximal Na/K pump rate
$G_\mathrm{glia,max}$	5 mM s−1	Maximal glial uptake strength of potassium
$\epsilon_\mathrm{K,max}$	0.25 s−1	Maximal potassium diffusion rate
$\epsilon_\mathrm{Na,max}$	0.1705 s−1	Maximal sodium diffusion rate
$\epsilon_\mathrm{Cl,max}$	0.2597 s−1	Maximal chloride diffusion rate
$\epsilon_\mathrm{O}$	0.17 s−1	Oxygen diffusion rate
K$^{+}_\mathrm{bath}$	6.25 mM	Normal bath potassium concentration
Na$^{+}_\mathrm{bath}$	147.4 mM	Normal bath sodium concentration
Cl$^{-}_\mathrm{bath}$	131.4 mM	Normal bath chloride concentration
α	5.3 g mol−1	Conversion factor
O$_{2 \mathrm{bath}}$	32 mg l−1	Normal bath oxygen concentration
U$_\mathrm{Kcc2}$	0.3 mM s−1	Maximal KCC2 cotransporter strength
U$_\mathrm{Nkcc1}$	0.1 mM s−1	Maximal NKCC1 cotransporter strength

The main premise of our temperature model is that cellular processes will accelerate when a moderate amount of heat is introduced. We used Q₁₀ factors to model these accelerated rates. The Q₁₀ factor of a process describes the change of the reaction rate when the temperature is changed by 10 ^∘C. For an arbitrary temperature difference, the Q₁₀ is calculated as:

$$\begin{align} Q_{10} = \left(\frac{y_{f}}{y_{0}}\right)^{10\ ^{\circ}\mathrm{C}/(\Delta T)} \end{align} \tag{ 9 }$$

where $\Delta T$ is the temperature change and y is the rate of a process. With a known Q₁₀, the rate at a new temperature T is given by:

$$\begin{align} y(T) = y_{0}Q_{10y}^{(T-T_{0})/(10\ ^{\circ}\mathrm{C})}. \end{align} \tag{ 10 }$$

Heat propagation is not instantaneous, and we assume that the temperature exponentially decays, given by:

$$\begin{align} T(t) = (T_{f}-T_{0})(1-e^{-t/\tau})+T_{0} \end{align} \tag{ 11 }$$

where τ is the time constant of the exponential decay. Combining equations (10) and (11) we get the time dependent rate equation:

$$\begin{align} y(t) = y_{0}Q_{10y}^{(T_{f}-T_{0})(1-e^{-t/\tau})/10}. \end{align} \tag{ 12 }$$

Using typical values reported in past literature (table 5), We estimated Q₁₀ model parameter values as reported in table 6. All full model computations were done in MATLAB (Mathworks).

Table 5. List of Q₁₀'s reported in literature.

Subject	Parameter	Q10
Giant squid axon [3]	Membrane conductivity	1.3
	Rate of change of gating variable	3
Rabbit sciatic nerve [86]	Time constant for h	1/3
	Conductivity of Na channels	1.7
Rat omohyoid muscle [87]	Time constant for n	1/2
Guinea pig hippocampal CA1 pyramidal neurons [88]	Inverse input resistance	1.72
Rat brain [89]	Immature NaKATPase activity	2.58
	Adult NaKATPase activity	3.05
Myelinated nerve fibres of Xenopus laevis [90]	$\alpha_{m},\beta{m}$	1.8, 1.7
	$\alpha_{h},\beta{h}$	2.8, 2.9
	$\alpha_{n},\beta{n}$	3.2, 2.8
	Sodium channel conductance	1.3
	Potassium channel conductance	1.2

Table 6.  Q₁₀'s used in the base model.

Parameter	Q10
$\frac{\mathrm{d}m}{\mathrm{d}t}$	1.75
$\frac{\mathrm{d}h}{\mathrm{d}t}$	2.85
$\frac{\mathrm{d}n}{\mathrm{d}t}$	3
$\rho_\mathrm{max}$	3
$G_\mathrm{Na}$, $G_\mathrm{K}$	1.7
$G_\mathrm{NaLeak}$, $G_\mathrm{ClLeak}$, $G_\mathrm{KLeak}$	1.1
$U_\mathrm{Nkcc1}$, $U_\mathrm{Kcc2}$	3

Bifurcations of the system's long-term behaviors were studied using numerical continuation software. We used MatCont [91], an open-source solver built on MATLAB. In order to improve reliability of this analysis, we used a reduced version of the model [43] where we replaced three ion concentrations with functions of other concentrations, as shown in equations (13). In addition, we also confined the cell volume to be constant.

$$\begin{align} [\mathrm{K}^{+}]_\mathrm{i}& = 140 + (18-[\mathrm{Na}^{+}]_\mathrm{i})-(6-[\mathrm{Cl}^{-}]_\mathrm{i})\nonumber \\[6pt] [\mathrm{Na}^{+}]_\mathrm{o}& = 144-\beta ([\mathrm{Na}^{+}]_\mathrm{i}-18)\nonumber\\[6pt] [\mathrm{Cl}^{-}]_\mathrm{o}& = 130-\beta ([\mathrm{Cl}^{-}]_\mathrm{i}-6) \end{align} \tag{ 13 }$$

Model parameteters are given in table 4

$$\begin{align} C\frac{\mathrm{d}V}{\mathrm{d}t}& = -I_\mathrm{Na}-I_\mathrm{K}-I_\mathrm{Cl}-I_\mathrm{pump}/\gamma \nonumber \\ \frac{\mathrm{NK}^+_\mathrm{o}}{\mathrm{d}t} & = \frac{1}{\tau}(\gamma\beta I_\mathrm{K} - 2\beta I_\mathrm{pump} -I_{k \mathrm{diff}} - I_\mathrm{glia}\nonumber\\ & \qquad\; - 2I_\mathrm{gliapump} +\beta I_\mathrm{KCC2} +\beta I_\mathrm{NKCC1})v_\mathrm{o} \nonumber \\ \frac{\mathrm{NK}^+_\mathrm{i}}{\mathrm{d}t} & = \frac{1}{\tau}(-\gamma I_\mathrm{K} + 2 I_\mathrm{pump} - I_\mathrm{KCC2} - I_\mathrm{NKCC1})v_\mathrm{i} \nonumber \\ \frac{\mathrm{NNa}^+_\mathrm{o}}{\mathrm{d}t} & = \frac{1}{\tau}(\gamma\beta I_\mathrm{Na}+3\beta I_\mathrm{pump}-I_\mathrm{Na diff}+\beta I_\mathrm{NKCC1})v_\mathrm{o} \nonumber \\ \frac{\mathrm{NNa}^+_\mathrm{i}}{\mathrm{d}t} & = \frac{1}{\tau}(-\gamma I_\mathrm{Na}-3 I_\mathrm{pump}- I_\mathrm{NKCC1})v_\mathrm{o}\nonumber\\ \frac{\mathrm{NCl}^-_\mathrm{o}}{\mathrm{d}t} & = \frac{1}{\tau}(-\gamma\beta I_\mathrm{Cl}-I_\mathrm{Cl diff}+\beta I_\mathrm{KCC2}+2\beta I_\mathrm{NKCC1})v_\mathrm{o} \nonumber \\ \frac{\mathrm{NCl}^-_\mathrm{i}}{\mathrm{d}t} & = \frac{1}{\tau}(\gamma I_{Cl}- I_\mathrm{KCC2}-2 I_\mathrm{NKCC1})v_\mathrm{o} \nonumber \\ \frac{[\mathrm{O}_2]}{\mathrm{d}t} & = -\alpha(I_\mathrm{pump}+I_\mathrm{gliapump})+\epsilon_O([\mathrm{O}_2]_\mathrm{bath}-[\mathrm{O}_2]) \nonumber \\ \gamma & = S/(F*v_\mathrm{i}) \nonumber \\ S & = 4\pi(3v_\mathrm{i}/4\pi)^{2/3} \nonumber \\ \tau& = 1000 \end{align} \tag{ 14 }$$

τ is the conversion factor from seconds to milliseconds, β is the ratio of intracellular volume to extracellular volume, S is the surface area of the cell, F is the Faraday constant, and γ is the conversion factor to convert concentration change (mM s⁻¹) into current density (µA cm⁻²).

The ionic currents and Nernst potentials are:

$$\begin{align} I_\mathrm{Na}& = G_\mathrm{Na}m^3h(V-E_\mathrm{Na})+G_\mathrm{NaL}(V-E_\mathrm{Na})\nonumber\\ I_\mathrm{K}& = G_\mathrm{K}n^4(V-E_\mathrm{K})+G_\mathrm{KL}(V-E_\mathrm{K})\nonumber\\ I_\mathrm{Cl}& = G_\mathrm{ClL}(V-E_\mathrm{Cl})\nonumber\\ E_\mathrm{K} & = 0.08617*T * \ln\left(\frac{[\mathrm{K}^+]_\mathrm{o}}{[\mathrm{K}^+]_\mathrm{i}}\right)\nonumber\\ E_\mathrm{Na} & = 0.08617*T * \ln\left(\frac{[\mathrm{Na}^+]_\mathrm{o}}{[\mathrm{Na}^+]_\mathrm{i}}\right)\nonumber\\ E_\mathrm{Cl}& = 0.08617*T * \ln\left(\frac{[\mathrm{Cl}^-]_\mathrm{i}}{[\mathrm{Cl}^-]_\mathrm{o}}\right) \end{align} \tag{ 15 }$$

where T is the temperature in K.

Rate of change of gating variables are given by:

$$\begin{align} \frac{\mathrm{d}m}{\mathrm{d}t}& = \alpha_{m}(1-m)-\beta_{m} \nonumber\\ \frac{\mathrm{d}h}{\mathrm{d}t}& = \alpha_{h}(1-h)-\beta_{h} \nonumber\\ \frac{\mathrm{d}n}{\mathrm{d}t}& = \alpha_{n}(1-n)-\beta_{n} \nonumber\\ \alpha_m & = \frac{0.32(54+V)}{1-e^{\frac{-(V+54)}{4}}} \nonumber\\ \beta_m & = \frac{0.28(V+27)}{e^{(V+27)/5}-1}\nonumber\\ \alpha_h & = 0.128e^{\frac{-(50+V)}{18}} \nonumber\\ \beta_h & = \frac{4}{1+e^{\frac{-(V+27)}{5}}} \nonumber\\ \alpha_n & = \frac{0.032(V+52)}{1-e^{\frac{-(V+52)}{5}}} \nonumber\\ \beta_n & = 0.5e^{\frac{-(V+57)}{40}}. \end{align} \tag{ 16 }$$

The pump currents and the glial currents are:

$$\begin{align} I_\mathrm{pump}& = \rho_{\max}*\frac{1}{1+\exp((20-[\mathrm{O}_2])/3)}\nonumber\\ & \quad *\frac{1}{1+\exp((25-[\mathrm{Na}^{+}]_\mathrm{i})/3)}\nonumber\\ & \quad *\frac{1}{1+\exp(3.5-[\mathrm{K}^{+}]_\mathrm{o})}\nonumber\\ I_\mathrm{gliapump}& = \frac{I_\mathrm{pump}}{3}\nonumber\\ I_\mathrm{glia}& = G_\mathrm{glia,max}*\frac{1}{1+\exp(18-[\mathrm{K}^{+}]_\mathrm{o})}\nonumber\\ & \quad *\frac{1}{1+\exp((2.5-[\mathrm{O}_2])/0.2)}. \end{align} \tag{ 17 }$$

Cell volume:

$$\begin{align} \hat{v}_\mathrm{i} & = v^{0}_\mathrm{i} * \left( 1.1029 - 0.1029 * \exp{\frac{\pi_\mathrm{o}-\pi_\mathrm{i}}{20}}\right)\nonumber\\ \pi_\mathrm{o} & = [\mathrm{Na}^+]_\mathrm{o} + [\mathrm{K}^+]_\mathrm{o} + [\mathrm{Cl}^-]_\mathrm{o} +[A^-]_\mathrm{o}\nonumber\\ \pi_\mathrm{i} & = [\mathrm{Na}^+]_\mathrm{i} + [\mathrm{K}^+]_\mathrm{i} + [\mathrm{Cl}^-]_\mathrm{i} +[A^-]_\mathrm{i}\nonumber\\ \frac{\mathrm{d}v_\mathrm{i}}{\mathrm{d}t}& = \frac{\hat{v}_\mathrm{i}-v_\mathrm{i}}{250} \end{align} \tag{ 18 }$$

where $[A^-]_\mathrm{i}$ and $[A^-]_\mathrm{o}$ are the concentrations of impermeable anions in the intracelluar and extracelluar space, respectively. $[A^-]_\mathrm{o}$ is assumed to be 18 mM while $[A^-]_\mathrm{i}$ is calculated at the initial condition such the osmotic pressure is at an equilibrium.

The non-electrogenic cotransporters NKCC1 and KCC2 are given by:

$$\begin{align} I_\mathrm{KCC2}& = L_\mathrm{KCC2}U_\mathrm{KCC2}\ln{\frac{[\mathrm{K}^+]_\mathrm{i}[\mathrm{Cl}^-]_\mathrm{i}}{[\mathrm{K}^+]_\mathrm{o}[\mathrm{Cl}^-]_\mathrm{o}}}\nonumber\\ I_\mathrm{NKCC1}& = U_\mathrm{NKCC2}((1-L_\mathrm{NKCC1})(f\,([\mathrm{K}^+]_\mathrm{o})-1)+1)\nonumber\\ & \quad \times \left(\ln{\frac{[\mathrm{K}^+]_\mathrm{i}[\mathrm{Cl}^-]_\mathrm{i}}{[\mathrm{K}^+]_\mathrm{o}[\mathrm{Cl}^-]_\mathrm{o}}}+\ln{\frac{[\mathrm{Na}^+]_\mathrm{i}[\mathrm{Cl}^-]_\mathrm{i}}{[\mathrm{Na}^+]_\mathrm{o}[\mathrm{Cl}^-]_\mathrm{o}}}\right)\nonumber\\ f\,([\mathrm{K}^+]_\mathrm{o})& = \frac{1}{1+e^{16-[\mathrm{K}^+]_\mathrm{o}}} \end{align} \tag{ 19 }$$

where $L_\mathrm{NKCC1}$ and $L_\mathrm{KCC2}$ are cotransporter levels.

Diffusion to the extracullar space is given by:

$$\begin{align} I_{x \mathrm{diff}} & = \epsilon_x([x]_\mathrm{o}-[x]_\mathrm{bath})\nonumber\\ \epsilon_x & = \epsilon_{x,\mathrm{max}}*\frac{1}{1+\exp((-20+\beta)/2)}\nonumber\\ & \quad * \frac{1}{1+\exp((2.5-[\mathrm{O}_2])/0.2)}\nonumber\\ x & = \mathrm{K,\,Na,\,Cl}. \end{align} \tag{ 20 }$$

The model was initialized by running the model from the initial conditions as shown in table A1 until the system reaches a steady-state.

Table A1. Coordinates used to initialize the model.

V	−73.59 mV
m	0.036
h	0.9992
n	0.0122
β	7
$[\mathrm{K}^+]_\mathrm{o}$	6.25 mM
$[\mathrm{K}^+]_\mathrm{i}$	140 mM
$[\mathrm{Na}^+]_\mathrm{o}$	147.4 mM
$[\mathrm{Na}^+]_\mathrm{i}$	18 mM
$[\mathrm{Cl}^-]_\mathrm{o}$	131.4 mM
$[\mathrm{Cl}^-]_\mathrm{i}$	6 mM
$[\mathrm{O}_2]$	29.3 mM
T	30 ∘C

An extended version of the $\dot{V}$ curve at the onset of the temperature step (figure 3(i)) is shown in figure A8. There exists three zeros at P1, P2, and P3 as shown, which are the equilibria of the system. Conventional methodolgy tells us that an equilibrium is stable if the slope of the first derivative is negative. We can see that P1 and P3 have negative slopes and P2 has a positive slope. While we can concude that P2 is an unstable equilibrium, this does not explain the full details for P1 and P3.

A true stable equilibrium will be stable against small perturbations to the gating variables. We tested the stability for small shifts in potential from $V_\mathrm{eq}$ to $V = V_\mathrm{eq}+\delta V$. The stability of the gating variables at the equilibria was determined by calculating the eigenvalues of the Jacobians under these assumptions. For stable equilibria, the eigenvalues must all have negative real parts. These calculations were done in Mathematica (Wolfram Research, Inc. Champaign, Illinois).

The system of equations for the voltage and the gating variables is the same as the full model's:

$$\begin{align} C\frac{\mathrm{d}V}{\mathrm{d}t}& = -I_\mathrm{Na}-I_\mathrm{K}-I_\mathrm{Cl}-I_\mathrm{pump}/\gamma \nonumber\\ \frac{\mathrm{d}m}{\mathrm{d}t}& = \alpha_{m}(1-m)-\beta_{m}\nonumber\\ \frac{\mathrm{d}h}{\mathrm{d}t}& = \alpha_{h}(1-h)-\beta_{h}\nonumber\\ \frac{\mathrm{d}n}{\mathrm{d}t}& = \alpha_{n}(1-n)-\beta_{n}. \end{align} \tag{ 21 }$$

Near some $V = V_\mathrm{i}$, the Jacobian is calculated as the matrix:

$$\begin{align*} \begin{bmatrix} \mathrm{d}v\mathrm{d}v_\mathrm{i} &\mathrm{d}v\mathrm{d}m_\mathrm{i} &\mathrm{d}v\mathrm{d}h_\mathrm{i}& \mathrm{d}v\mathrm{d}n_\mathrm{i}\\ \mathrm{d}m\mathrm{d}v_\mathrm{i} & \mathrm{d}m\mathrm{d}m_\mathrm{i} & 0 & 0\\ \mathrm{d}h\mathrm{d}v_\mathrm{i} &0 &\mathrm{d}h\mathrm{d}h_\mathrm{i}& 0\\ \mathrm{d}n\mathrm{d}v_\mathrm{i} & 0 & 0 & \mathrm{d}n\mathrm{d}n_\mathrm{i}\\ \end{bmatrix} \end{align*}$$

where the elements are given by:

$$\begin{align} \mathrm{d}v\mathrm{d}v_\mathrm{i} & = \frac{\partial}{\partial V} \left. \left( \left. \frac{\partial V}{\partial t} \right \rvert_{m = m_{\infty}(V_\mathrm{i}),h = h_{\infty}(V_\mathrm{i}),n = n_{\infty}(V_\mathrm{i})} \right) \right \rvert_{V = V_\mathrm{i}}\nonumber\\ \mathrm{d}v\mathrm{d}m_\mathrm{i} & = \frac{\partial}{\partial m} \left. \left( \left. \frac{\partial V}{\partial t} \right \rvert_{h = h_{\infty}(V_\mathrm{i}),n = n_{\infty}(V_\mathrm{i})} \right) \right \rvert_{V = V_\mathrm{i},m = m_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}v\mathrm{d}h_\mathrm{i} & = \frac{\partial}{\partial n} \left. \left( \left. \frac{\partial V}{\partial t} \right \rvert_{m = m_{\infty}(V_\mathrm{i}),n = n_{\infty}(V_\mathrm{i})} \right) \right \rvert_{V = V_\mathrm{i},h = h_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}v\mathrm{d}n_\mathrm{i} & = \frac{\partial}{\partial h} \left. \left( \left. \frac{\partial V}{\partial t} \right \rvert_{m = m_{\infty}(V_\mathrm{i}),h = h_{\infty}(V_\mathrm{i})} \right) \right \rvert_{V = V_\mathrm{i},n = n_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}m\mathrm{d}v_\mathrm{i} & = \frac{\partial}{\partial V} \left. \left( \frac{\partial m}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},m = m_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}h\mathrm{d}v_\mathrm{i} & = \frac{\partial}{\partial V} \left. \left( \frac{\partial h}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},h = h_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}n\mathrm{d}v_\mathrm{i} & = \frac{\partial}{\partial V} \left. \left( \frac{\partial n}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},n = n_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}m\mathrm{d}m_\mathrm{i} & = \frac{\partial}{\partial m} \left. \left( \frac{\partial m}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},m = m_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}h\mathrm{d}h_\mathrm{i} & = \frac{\partial}{\partial h} \left. \left( \frac{\partial h}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},h = h_{\infty}(V_\mathrm{i})}\nonumber\\ \mathrm{d}n\mathrm{d}n_\mathrm{i} & = \frac{\partial}{\partial n} \left. \left( \frac{\partial n}{\partial t} \right) \right \rvert_{V = V_\mathrm{i},n = n_{\infty}(V_\mathrm{i})}. \end{align} \tag{ 22 }$$

For point 1, all eigenvalues were real and negative, making this equilibrium stable. Point 3 also had all real eigenvalues, but only 2 of them were negative, thus making it unstable to small changes in the gating variables.

Analysis of the field map around these points confirms that point 1 is the only true stable point. Near these equilibria, $\mathrm{d}m/\mathrm{d}t \gg \mathrm{d}n/\mathrm{d}t \gg \mathrm{d}h/\mathrm{d}t$ (figure A9). Thus, we assume that for each gating variable, faster gating variables are at their steady-state values for the corresponding voltage ($x_{\infty} (V)$), and slower gating variables are at their steady-state values for the equilibrium voltage ($x_{\infty} (V_\mathrm{eq})$).