Hypersynchronous (HYP) seizure onset is one of the frequently observed seizure-onset patterns in temporal lobe epileptic animals and patients, often accompanied by hippocampal sclerosis. However, the exact mechanisms and ion dynamics of the transition to HYP seizures remain unclear. Transcranial magneto-acoustic stimulation (TMAS) has recently been proposed as a novel non-invasive brain therapy method to modulate neurological disorders. Therefore, we propose a biophysical computational hippocampal network model to explore the evolution of HYP seizure caused by changes in crucial physiological parameters and design an effective TMAS strategy to modulate HYP seizure onset. We find that the cooperative effects of abnormal glial uptake strength of potassium and excessive bath potassium concentration could produce multiple discharge patterns and result in transitions from the normal state to the HYP seizure state and ultimately to the depolarization block state. Moreover, we find that the pyramidal neuron and the PV+ interneuron in HYP seizure-onset state exhibit saddle-node-on-invariant-circle/saddle homoclinic (SH) and saddle-node/SH at onset/offset bifurcation pairs, respectively. Furthermore, the response of neuronal activities to TMAS of different ultrasonic waveforms revealed that lower sine wave stimulation can increase the latency of HYP seizures and even completely suppress seizures. More importantly, we propose an ultrasonic parameter area that not only effectively regulates epileptic rhythms but also is within the safety limits of ultrasound neuromodulation therapy. Our results may offer a more comprehensive understanding of the mechanisms of HYP seizure and provide a theoretical basis for the application of TMAS in treating specific types of seizures.

HYP seizure onset is characterized as low-frequency, high-amplitude periodic spikes and is thought to result from an imbalance between inhibition and excitation. As a novel neuromodulation, TMAS has attracted widespread attention for its non-invasive regulation of neural activity without tissue damage. However, due to neurons’ intrinsic characteristics and interactions, different stimulation parameters may lead to transitions in neuronal states. Computational models can reproduce phenomena observed from physiological experiments, reveal underlying seizure mechanisms and neuromodulation from the perspective of dynamics, and guide future treatment for epileptic disorders. We complement our network model with variables for the ion concentration to reveal the dynamical mechanisms of the evolution of HYP seizure driven by excess extracellular potassium. Furthermore, we design an effective TMAS modulation strategy to control HYP seizures in the human hippocampus.

## I. INTRODUCTION

Temporal lobe epilepsy in humans is caused by abnormal neural discharges in the focus lesion of the brain, which is often closely associated with the hippocampus. According to electroencephalographic recordings in temporal lobe epileptic patients and animal models, hypersynchronous (HYP) and low-voltage fast (LVF) are identified as the most frequently observed focal seizure-onset patterns, and the LVF onset pattern often follows HYP onset.^{1–3} High-amplitude periodic spikes that occur at a frequency below 2 Hz are considered electrophysiological characteristics of HYP onset seizure.^{4} Previous studies and experiments have shown that HYP onset seizure is caused by progressive exhaustion of inhibition and unrestrained enhancement of excitation.^{5–7} In recent years, it has been documented that the HYP onset pattern can be triggered by optogenetic activation of pyramidal neurons in the *vitro* 4-aminopyridine model and is mostly accompanied by fast ripples (250–500 Hz).^{8,9} With the increasing research on neurological diseases, the role of ions in temporal lobe epilepsy is receiving more attention. From the perspective of ion dynamics, there are several possibilities for its contribution to epilepsy: for instance, extracellular potassium accumulation due to blockade of Kv1 potassium channels^{10} or impaired astrocytic potassium regulation system;^{11} low magnesium;^{5} excess intracellular chloride due to abnormal activities of $ N a +/ K +/2 C l \u2212$ cotransporter-1 (NKCC1) and $ K +/ C l \u2212$ cotransporter-2 (KCC2) cotransporter.^{12}

There are numerous mathematical models proposed to study the effect of ion dynamics on the transitions from regular neuronal activity to epilepsy. Ullah presented a model containing the following compartments: neurons, glia, and extracellular space to investigate the effects of high potassium, glial cell buffering,^{13} and cell swell^{14} on epilepsy. Furthermore, Wei *et al.*^{15} used a similar model to reproduce LVF seizure-like events caused primarily by high potassium due to hypoxia and validated the *vitro* experiments in the Ammon's horn 1 (CA1) region of rat hippocampal slices.^{16,17} In addition, Gentiletti constructed a small-scale biophysical network to explore neuronal state transitions between the various LVF seizure phases and corresponding changes in Na, K, and Cl ion concentrations.^{18} Yao used a single-neuron model to study the dynamic mechanisms of seizure and spreading depression based on oxygen and ion concentrations.^{19} A minimal model called “Epileptor-2” was proposed to indicate that ictal discharges were triggered by potassium accumulation and eliminated by sodium–potassium pump activation.^{20} However, to the best of our knowledge, in contrast to numerous computational models of LVF onset seizure, only a few were proposed to study mechanisms of HYP seizure onset at the ion concentration level. In particular, Buchin *et al.*^{21} presented a hippocampus subiculum network model to study the transition from gamma oscillations to HYP seizure due to the reduced efficacy of the KCC2 cotransporter, but this research mainly focused on transitions between normal and epileptic state rather than on the complete dynamic state transitions.

The existing treatment methods for epilepsy mainly include medication, surgical treatment, and neuromodulation. However, medication lacks brain target specificity, and surgical treatment by implanting electrodes can cause irreversible damage to the brain. Compared with traditional brain stimulation such as deep brain stimulation, optogenetic stimulation, and transcranial magnetic stimulation, transcranial magneto-acoustic stimulation (TMAS) has the advantage of non-invasiveness, high spatial resolution (<3 mm), and high penetration depth.^{22,23} Therefore, TMAS is gradually playing an important role in treating neurological diseases such as Alzheimer's disease^{24} and Parkinson's disease.^{25} As a novel non-invasive neuromodulation, TMAS can generate a stimulating electric current in nerve tissues by Lorentz forces on moving ions produced by the ultrasonic wave in the presence of a magnetic field.^{26} Previous physiological experiments have shown that low-intensity pulsed ultrasound stimulation can effectively regulate neuronal potassium currents^{27} and show great potential in the suppression and treatment of epilepsy.^{28,29} Moreover, Yuan *et al.*^{30,31} demonstrated that TMAS could affect neuronal firing rhythms and ion concentrations, which provided a theoretical basis for TMAS in the modulation of HYP seizures. However, the crucial ultrasound stimulation parameters and therapeutic effects of TMAS on epilepsy need further research.

Thus, in this paper, we construct a biophysical network model containing pyramidal neurons, a parvalbumin-positive (PV+) interneuron, and astrocytes to investigate the inherent ionic mechanisms of the complete evolution of HYP seizure and the dynamical properties of bifurcation types at seizure onset/offset. Further, we explore precise control of the HYP epileptic hippocampal network by designing an effective TMAS regulation strategy. Our results may provide a theoretical basis and guidance for diagnosing and treating HYP seizures.

## II. METHODS AND MODEL

### A. Hippocampal network model

Our biophysical computational hippocampal Ammon's horn 3 (CA3) network model comprises four pyramidal neurons and a PV+ interneuron. Each neuron is embedded with extracellular space and an astrocyte. The extracellular space of each neuron is enclosed in a bath solution that represents blood vessels *in vivo*. As shown in Fig. 1(a), pyramidal neurons receive excitation from other pyramidal neurons via glutamatergic receptors and inhibition from a PV+ interneuron via gamma-aminobutyric acid type A $( GAB A A)$ receptor, whereas the PV+ interneuron only receives glutamate from pyramidal neurons. It is worth mentioning that our network connectivity is all-to-all, so there are synaptic connections between the pyramidal neurons located on the diagonal of Fig. 1(a) (PY1–PY3 and PY2–PY4). For instance, the four lines on the main diagonal that involve PY3 in Fig. 1(a) represent the synaptic connections of PY3 to PY1, PY3 to PV+, PV+ to PY3, and PY1 to PY3.

^{32}PV+ basket neurons are thought to account for 60% of all PV-expressing interneurons, which are the most significant number of neurons in all PV-expressing interneurons.

^{33}Thus, we used PV+ basket interneurons modeled by Wang and Buzsáki,

^{34}as shown in Eq. (2). More details on neuron modeling are available in Refs. 31 and 33,

*C*,

*V*, and

*t*denote capacitance density, neuronal membrane potential, and time, respectively. The units of the above variables are $\mu F/ c m 2$, mV, and ms. The pyramidal neuron contains sodium voltage-gated current $ I Na$, persistent sodium voltage-gated current $ I Nap$, potassium voltage-gated current $ I K$, calcium-activated potassium or after-hyperpolarization (AHP) current $ I AHP$, leak current $ I L$, synaptic current $ I syn$, and sodium–potassium pump current $ J pump/\gamma $, whereas the PV+ interneuron does not contain $ I Nap$ and $ I AHP$. $\gamma $ represents the conversion factor from the current density into the ion concentration rate of change. For the pyramidal neuron, we set $ \gamma PY=0.0445 mM/( C\u22c5 cm)$. For the PV+ interneuron, we set $ \gamma PV +=0.0286 mM/( C\u22c5 cm)$.

^{35}The voltage-gated ion channel equations are shown in Eq. (3),

^{36}

^{,}

*x*= Na, K, and Cl) can be calculated from the dynamic intracellular $( [ x ] i)$ and extracellular $( [ x ] o)$ ionic concentrations by the Nernst equation. In particular, GABA receptors are approximately four times more permeable to chloride than bicarbonate.

^{37}Thus, the effect of chloride and bicarbonate ions on the GABAergic reversal potential $ E GABA$ in our model is set to four to one,

*x*= GABA, AMPA, or NMDA and the synaptic current is from neuron

*i*to neuron

*j*. $ s x , i$ (

*x*= GABA, AMPA, and NMDA) represents the fraction of synaptic gating variable related to the presynaptic neuron

*i*. $ [ Mg ] o$ and $B( V j)$ denote extracellular magnesium concentration and postsynaptic membrane potential dependence fraction on magnesium, respectively. We set $ [ Mg ] o=2 mM$. $ \tau R$ and $ \tau D$ represent the synaptic rise and decay time constant, respectively. The AMPAergic reversal potential $ E AMPA$ and the NMDAergic reversal potential $ E NMDA$ are set to zero. For an AMPA-type synapse, we set $ \tau R=0.1$, $ \tau D=3$, $ g PY \u2192 PV +=0.01 mS/ c m 2$, and $ g PY \u2192 PY=0.04 mS/ c m 2$. For an NMDA-type synapse, we set $ \tau R=13$, $ \tau D=150$, $ g PY \u2192 PV +=0.3 mS/ c m 2$, and $ g PY \u2192 PY=3 mS/ c m 2$. For a GABA-type synapse, we set $ \tau R=0.3$, $ \tau D=40$, and $ g PV + \u2192 PY=0.15 mS/ c m 2$.

*et al.*shown in Eq. (6),

^{38}

### B. Ionic dynamics

We complemented the Hodgkin–Huxley formulation model with the ion concentration using differential equations based on previous studies.^{39} The specific details are shown in Fig. 1(b). The ionic dynamics of $ N a +$, $ K +$, $ C l \u2212$, $ C a 2 +$, and $ HCO 3 \u2212$ are regulated by ion channels, cotransporters (KCC2 and NKCC1), and synaptic channels. The intracellular and extracellular ion concentrations of each ion type are continuously updated as slow variables.

^{39}

^{,}

^{40}The mean distance between two neurons $dx=50\mu m$ is measured in the hippocampus.

^{12}$\epsilon $ and $ K bath$ denote the maximal potassium diffusion rate and potassium concentration from a bath solution, respectively. We set parameter $ \rho pump=1.25 mM/ s$, $ U nkcc 1=0.1 mM/ s$, and $ U kcc 2=0.3 mM/ s$.

### C. Theoretical basis for TMAS

^{31}Schematic diagram of TMAS is shown in Fig. 2(a). Therefore, based on the previous works,

^{30,41}the TMAS-induced electric current density $ I TMAS$ can be added directly to Hodgkin–Huxley models expressed as

*f*represent nerve tissue conductance, tissue density, ultrasound speed, and ultrasound frequency, respectively. Magnetic field intensity

*B*and ultrasound intensity $ U I$ will be further analyzed as crucial modulation parameters. The bias parameter

*b*determines the ultrasonic waveforms, which we describe in detail in Fig. 2(b). Here, we set $\sigma =0.5 Simens/ m$, $\rho =1120 Kg/ m 3$, $ c 0=1540 m/ s$, and $f=0.2 MHz$.

*MI*) indicates the potential biomechanical danger due to ultrasound cavitation. The thermal index (

*RI*) reflects the temperature rise in neural tissue due to the absorption of ultrasound energy. With reference to standard regulation of the United States Food and Drug Administration (FDA),

^{43}combined with our model, the MI and RI formulas are simplified as shown in Eqs. (12) and (13), respectively,

*S*is the ultrasound spot area in cm

^{2}, $ C MI=1 MPa\u22c5 MH z \u2212 1 / 2$, and $ C TIS=210 mW\u22c5 MHz$.

All simulations are carried out in MATLAB, version R2021b (https://ww2.mathworks.cn). Differential equations are solved with the fourth-order Runge–Kutta method. Taking the high-frequency characteristics of ultrasound into consideration, we set a fixed time step $t=0.0005 ms$. Taking the time step and ultrasound frequency into Eq. (10), we obtain the sine function with a time step $2\pi ft=0.2\pi $. Ten time steps are used to simulate a complete sine cycle, which ensures the precision of TMAS. The Parallel Computing Toolbox in MATLAB is used to speed up computation when we solve equations in parameter space. In addition, the initial conditions for membrane potentials and ion concentrations are determined based on the normal distribution with the following means: $ V PY=\u221267 mV$, $ V PV +=\u221265 mV$, $ [ K ] o=4 mM$, $ [ Na ] i=18 mM$, $ [ Cl ] i=6 mM$, and variance at 5% of the corresponding mean.

## III. RESULTS

### A. Dynamic transition to HYP seizures due to abnormal uptake of extracellular potassium

In this section, we focus on the complete state transitions and the corresponding ionic concentration dynamics due to changes in bath potassium concentration $( K bath)$ and the glial uptake strength of potassium $( G glia)$ parameters which may be associated with HYP seizures observed in hippocampal slices in high $ K +$ medium.^{44} Alterations in physiologic parameters $ K bath$ and $ G glia$ cause rapid changes in $ [ K ] o$ within pyramidal neurons and PV+ interneuron, accompanied by subsequent changes in the other ion concentrations $( [ Na ] i , [ Cl ] i , [ Ca ] i)$, and ultimately induction of state transitions. As shown in Fig. 3, two-parameter space diagrams for parameters $ K bath$ and $ G glia$ are elaborated to study the quantitative effects on the network's state transitions. There are six states distinguished by discharge amplitude and interspike interval (ISI) including the resting state (1), fast-spiking state (2), synchronous bursting state (3), HYP seizure-onset state (4), sustained ictal activity state (5), and sustained depolarization block state (6). We can see that elevating $ K bath$ can gradually cause neurons to transition from normal discharges (resting state or fast-spiking state) to HYP seizures (synchronous bursting state, HYP seizure-onset state or sustained ictal activity) and eventually to the depolarization block state. However, increasing $ G glia$ improves the glial uptake strength of potassium and effectively reduces $ [ K ] o$, converting HYP seizure onset into synchronous bursting discharges or fast-spiking state.

To further understand the role of $ K bath$ in the neural discharge patterns and ionic concentration dynamics in the network model, we gradually increase $ K bath$ and fix *G*_{glia} = 10 mM/s. As a result of GABA inhibitory signaling from the interneuron, all four pyramidal neurons have essentially the same discharge pattern after entering the steady state. Therefore, we only present one of the pyramidal neurons to illustrate the differences between the two types of neurons in different states. Figure 4 shows the time series of the membrane potentials and ion concentrations in two classes of neurons vs different parameters $ K bath$. As shown in Fig. 4(a1), all neurons remain at a resting state and ion concentrations remain stable when *K*_{bath} = 4 mM/s. As $ K bath$ increases, slight increases in $ [ K ] o$ make the neurons more excitable, resulting in the network transition to a fast-spiking state [Fig. 4(a2)]. According to the first and third panels in Fig. 4(a2), it can be clearly seen that the PV+ interneuron discharges faster than the pyramidal neuron. When $ K bath$ further increases, the pyramidal neuron and PV+ interneuron appear in a synchronized bursting state [Fig. 4(a3)], and their ion concentrations begin to fluctuate dramatically within a certain range. When $ K bath$ approximately approaches 16 mM, the pyramidal neuron remains bursting, while the PV+ interneuron partially enters a phase of depolarization block [Fig. 4(a4)]. This phenomenon is defined as the HYP seizure-onset state, which is similar to physiological experiments reported by Zhang *et al.*^{6} and Cammarota *et al.*^{45} The phase of depolarization block in PV+ interneuron is mediated by the inactivation of sodium channels due to extracellular potassium accumulation, which leads to a decrease in inhibitory GABA currents and more excitation of pyramidal neurons. Further increasing $ K bath$, the neural network is considered to enter a sustained ictal activity [Fig. 4(a5)]. In this state, the membrane potentials and ion concentrations of the pyramidal neuron enter a state of high-frequency slight oscillation. However, the interneuron remains in the sustained depolarization block and its extracellular potassium ion concentration remains constant at a high value. Eventually, all neurons enter a complete depolarization block state when $ K bath$ reaches a higher level [Fig. 4(a6)].

We further analyze the HYP seizure mechanism from the perspective of nonlinear dynamics combined with ion concentration dynamics. Complex oscillations such as bursting are triggered by fast/slow systems activity, where the slow subsystem drives the fast subsystem to transition between different states.^{46} In our model, the neuronal membrane potential timescale is milliseconds as a variable in the fast subsystem, while the ion concentration timescale is seconds as a variable in the slow subsystem. Figure 5 illustrates the evolution of the three-dimensional phase space of the two types of neurons in different firing patterns. From a dynamical perspective, neuronal resting and fast-spiking states are monostable. They correspond to a stable fixed point [Fig. 5(a1)] and a limit cycle [Fig. 5(a2)] in the phase space, respectively. The transition from the resting state to the fast-spiking state involves a saddle-node-on-invariant-circle (SNIC) bifurcation. After entering HYP seizures, the system switches between the quiescent and oscillatory behaviors. Jirsa *et al.* proposed a dynamical taxonomy of seizure-like events that systematically classified epileptic seizures into 16 different onset/offset bifurcation pairs.^{47} Saggio *et al.* provided a dynamical framework that can easily identify the features of oscillation onset/offset bifurcation pair by neuronal discharge amplitude, ISI, baseline jump, and flow changes in the phase space.^{48,49} As shown in Fig. 5(a3), the phase space of the pyramidal neuron starts at a stable fixed point corresponding to a quiescent state. The reduction of slow-variable $ [ K ] o$ drive trajectory to slide along the curve with a black arrow, where the black arrow indicates the direction of movement. When $ [ K ] o$ approximately approaches 13 mM, the potassium gate variable is activated and the trajectory transits to the stable limit cycle corresponding to the oscillatory state by crossing an SNIC onset bifurcation. Then, $ [ K ] o$ begins to gradually increase and the trajectory displays the large amplitude oscillations forming a tube-like structure. As $ [ K ] o$ approximately approaches 15 mM, the system reaches an offset bifurcation and the oscillation ends. At this time, the trajectory jumps back to a stable fixed point by crossing an SH offset bifurcation. According to the dynamical taxonomy of Izhikevich, this would represent a “circle/homoclinic” bursting.^{46} However, the phase space of the PV+ interneuron seems to be different. The trajectory starts from a stable point and undergoes small oscillations triggered by the excitation of pyramidal neurons. Then, the trajectory jumps up to the stable limit cycle by crossing an SN onset bifurcation. This leads to a baseline jump in the amplitude of PV+ interneuron from a quiescent state to an oscillatory state [third panel in Fig. 4(a3)]. After crossing the onset bifurcation, the trajectory exists in the form of limit cycles. Eventually, the trajectory jumps back to a stable fixed point by crossing an SH offset bifurcation. This leads to a sudden jump back to the baseline from an oscillatory state to a quiescent state [third panel in Fig. 4(a3)]. Similarly, the pyramidal neuron and PV+ interneuron in HYP seizure-onset state show SNIC/SH and SN/SH bifurcations pair at seizure onset/offset, respectively [Fig. 5(a4)]. It is worth noting that the trajectory of PV+ interneuron in HYP seizure-onset state shrinks from the limit cycle to a point corresponding to depolarization block via a supercritical Hopf (SupH) offset when $ [ K ] o$ approximately approaches 15.2 mM. Then, the neuronal amplitude increases again and the trajectory transit from the point to the limit cycle via a supercritical Hopf (SupH) offset. However, the final bifurcation offset is SH. Hence, we still classify this bifurcation pair type as an SN/SH. When $ [ K ] o$ increases beyond a certain threshold, the neurons transition to a depolarizing block state, causing the system to re-enter monostability via SupH bifurcation. It is just like the phase space of the PV+ interneuron shown in Fig. 5(a5) and two types of neurons in Fig. 5(a6).

### B. The therapeutic effect of TMAS on the HYP seizure onset under different ultrasonic parameters

In this section, we investigate the therapeutic effect of TMAS on HYP seizure onset where different ultrasonic waveforms are chosen to target the pyramid neurons and PV+ interneuron. After the system is in a homeostatic HYP seizure onset for 7 s, continuous TMAS currents are simultaneously applied to both types of neurons. Figures 6(a1) and 6(a2) show the curve change of the pyramidal neurons and PV+ interneuron firing rates (FR) under different ultrasound waveforms with ultrasound intensity on HYP seizure onset, respectively. We find that the pyramidal neuron and PV+ interneuron firing rates gradually decrease with the increase of $ U I$ under lower sine wave stimulation, even reaching zero when ultrasound intensity exceeds 0.5 W/cm^{2}. However, pyramidal neuron and PV+ interneuron firing rates gradually increase under upper sine wave stimulation. Interestingly, two types of neuronal firing rates are basically unchanged under sine wave stimulation and remain similar to those in the absence of stimulation. In Figs. 6(b1)–6(b3), we present the time series of the membrane potentials and ion concentrations in two classes of neurons under the upper lower wave, sine wave, and upper sine wave stimulation when the parameters $B=0.9 T$ and $ U I=0.7 W/ c m 2$. The lower sine wave can inhibit neuron discharges and lead to the network transition from HYP seizure-onset state to rest state. In contrast, the upper sine wave can promote neuron discharges and increase the duration of HYP seizure. The sine wave seems to have no effect on neuronal activity, possibly because the promoting and inhibiting effects of the positive and negative phases of the sine wave cancel each other at the microsecond level, and it is ineffective in changing membrane potentials (at the millisecond level) and ion concentrations (at the second level) on a much larger time scale.

After exploring the previous results, we choose the continuous lower sine wave ultrasound as the desired TMAS waveform. Below, we examine the threshold of different stimulation parameters for the HYP seizure onset. Compared with the neuronal firing rate which mainly indicates the overall activity level and responsiveness of neurons, ISI can accurately provide information about changes in individual neurons’ discharge frequency and bifurcations of discharge patterns from a dynamical perspective. Therefore, we calculate the ISI of neurons as a metric to study the effect of TMAS on neuronal spike frequency at different magnetic field intensities. It is observed from Fig. 7 that the larger value of the ISI (>1 s) represents the interval between two seizures, while the smaller value represents the time interval between two neuronal spikes in a seizure. Therefore, we hypothesize that the maximal ISI of pyramidal neurons after TMAS is an indicator of epileptic latency. As shown in Fig. 7, we find the maximum ISI of the pyramidal neuron and interneuron gradually increases before reaching the bifurcation point where $B=0.85 T$. It is worth mentioning that the maximum ISI of pyramidal neurons exceeds five seconds when $0.6 T\u2264B\u22640.8 T$. After crossing the bifurcation point, the ISI of the pyramidal neuron and interneuron suddenly drops to zero. We quantitatively identified a bifurcation point for state transitions by varying the ultrasound parameter *B*. These results indicate that lower sine wave TMAS can increase the latency of HYP seizures and eventually result in the transition to the resting state at sufficient intensity.

Next, we investigate the state transitions and the corresponding ionic concentration dynamics due to changes in crucial ultrasound parameters including *B* and $ U I$ under lower sine wave stimulation. The dependence of the control effect on *B* and $ U I$ is shown in Fig. 8, distinguished by neuronal discharge amplitude and ISI. Here, the red parameter area indicates epileptiform discharges and is considered as an invalid control. It can be observed in Fig. 8(b1) that the neurons still undergo HYP seizures due to insufficient stimulation when $ U I=0.1 W/ c m 2$ and $B=0.9 T$. As $ U I$ and *B* increase together, the latency of HYP seizures gradually exceeds 5 s, and, therefore, the yellow parameter area is classified as a partially invalid control. Compared with those in Fig. 8(b1), the latency of HYP seizure in Fig. 8(b2) significantly increases and HYP seizure duration in Fig. 8(b2) significantly decreases at the same regulatory time. When *B* and $ U I$ are further increased, the blue parameter area shows that neurons after TMAS enter the resting state, which is considered as a valid control. As shown in Fig. 8(b3), the membrane potentials of pyramidal neuron and interneuron were essentially stabilized at −67 and −70 mV, respectively, similar to the membrane potentials of neurons at rest rather than depolarization block. Meanwhile, we find decay of $ [ Na ] i$ and $ [ K ] o$ in two classes of neurons, which may result from activation of the $ N a +/ K +$ pump.

## IV. DISCUSSION AND CONCLUSION

In this paper, we propose a biophysical network to study state transitions from the normal state to the HYP seizure state (and ultimately to the depolarization block state) and explore the oscillatory dynamics of ion concentrations in different states, which is the first time in the model network to reproduce the complete evolution of HYP seizure driven by excess extracellular potassium. The results imply that HYP seizures are triggered by a dynamic decrease in inhibitory $ GAB A A$ signaling due to excess bath potassium concentration and abnormal glial uptake strength of potassium. Moreover, we find that the pyramidal neuron and PV+ interneuron in HYP seizure-onset state exhibit SNIC/SH and SN/SH bifurcation pairs at seizure onset/offset, respectively. On this basis, we study the therapeutic effects of TMAS in modulating HYP seizure onset under different ultrasonic waveforms and eventually find that negative waves can effectively increase the latency of HYP seizures and even completely suppress epilepsy by the metric of firing rate. Further, we focus on the modulation of HYP seizure-onset thresholds by lower sine wave TMAS and find that the greater the magnetic field strength and ultrasonic intensity, the better the control effect. More importantly, we show the parameter area that efficiently modulates HYP seizure onset.

Next, we discuss two ionic factors that may contribute to HYP seizure. First, the PV+ interneuron gradually enters a phase of depolarization block as $ [ K ] o$ exceeds 16 mM, which leads to more excitation of pyramidal neurons due to the lack of inhibitory GABA currents. When the interneuron enters the sustained depolarization block, the system converts to the sustained ictal activity state because the pyramidal neurons only receive excitatory glutamate signals from other pyramidal neurons rather than inhibitory GABA signals from the PV+ interneuron. *In vitro* studies in the isolated rat hippocampus in 4-AP or low-Mg medium have also shown that pyramidal neurons are recruited to ictal discharges when GABAergic interneurons undergo a depolarizing block phase.^{16,45} Interestingly, in our model, the interneuron is more likely to enter a depolarizing block than pyramidal neurons, which is also consistent with physiological experiments.^{50} For this phenomenon, we are inclined to hypothesize that the AHP current in pyramidal neurons prevents the depolarization block. Second, HYP seizures may be related to intracellular chlorine accumulation. Work in temporal lobe slices indicated that cotransporter NKCC1 activity is upregulated and cotransporter KCC2 activity is downregulated in epileptic patients with hippocampal sclerosis.^{12,51} In our study, it can be seen in Fig. 4 that chloride ions gradually influx during seizures in response to cotransporter NKCC1 and KCC2 activities. Then, according to Eq. (4), we can find a gradual positive shift in the GABA_{A} reversal potential, which is closely linked to postsynaptic GABA_{A} receptor activation in pyramidal neurons, resulting in dynamically weakening inhibitory $ GAB A A$ signaling. In addition, Weiss suggested that HYP onset may be caused by depolarizing inhibitory postsynaptic currents due to chloride ion dysregulation.^{52} Furthermore, intracellular recordings in rat perirhinal cortices during continuous 4-AP application indicated that HYP ictal onset may be triggered by diminished GABA_{A} signaling due to extracellular potassium accumulation.^{10} The above two studies validate our proposed mechanism of HYP seizure onset.

Numerous biophysical models using potassium ions as a slow variable in dynamic analyses have shown that the single pyramidal neuron transitions to seizure through a SNIC bifurcation and transitions to depolarization block through a Hopf bifurcation.^{53,54} However, many phenomenological models with few variables and parameters are characterized by SN/SH bifurcation pairs at the seizure onset/offset.^{47,55} In our model, the pyramidal neuron and the PV+ interneuron in the HYP seizure-onset state exhibit SNIC/SH and SN/SH bifurcation pairs at seizure onset/offset, respectively. Meanwhile, the behaviors of the pyramidal neuron and PV+ interneuron are similar to seizure-like events dynamics observed in a simplified single-neuron model,^{56} corresponding to paths “c6” and “c2” in existing dynamic frameworks, respectively.^{48} We emphasized the bifurcation of onset and offset pairs at specific stages of HYP seizure based on invariant dynamics characterizations, which provide new insight into epileptic mechanisms and treatment. However, due to random noise and ambiguous data, the classifications of epileptic bifurcations in clinical conditions have been a challenge. Therefore, effective identification and quantification of dynamic features may require more work in the future.

As a novel neuromodulation, the safety of TMAS must be taken into consideration. Previous research has shown that high-intensity focused ultrasound causes irreversible nerve tissue necrosis due to thermal injury caused by prolonged exposure, whereas lower-intensity focused ultrasound can activate or inhibit bioelectric activity without any concomitant brain damage.^{57,58} Therefore, three medical ultrasound metrics including spatial-peak temporal-average intensity $( I spta)$, *MI*, and *RI* are calculated to estimate the possible biological effects of ultrasound. It is worth noting that the US FDA allows diagnostic systems with an upper limit of 0.72 W/cm^{2} in $ I spta$, 1.9 in *MI*, and 6 in *RI*, which is also used as a reference for ultrasound neuromodulation parameters.^{59} Hence, in our model, we set the maximum $ U I$ to 0.7 W/cm^{2} to minimize possible thermal effects and comply with diagnostic regulations of $ I spta$. We can calculate $ p L=0.155 MPa$, $S=7.065\xd7 10 \u2212 2 c m 2$ with the ultrasonic parameters $ U I=0.7 W/ c m 2$, $f=0.2 MHz$, and the diameter of the ultrasonic spot equals to $3 mm$. Taking the above results into Eqs. (11) and (12), we obtain $MI=0.347$ and $RI=0.047$, which are well below the FDA recommended limit that can cause thermal or mechanical damage to nerve tissue. Meanwhile, for the sake of safety, our magnetic field intensity $( B \u2264 1 T)$ is also within the parameter range required for diagnosis. In addition, ultrasound frequency is also a safety indicator of great concern. The frequency range of diagnostic medical ultrasound is between 1 and 15 MHz, while the frequency range of therapeutic ultrasound is in kHz.^{60,61} Low-frequency ultrasound (*f* < 0.65 MHz) can effectively reduce ultrasonic attenuation in the scalp and skull and is widely used to stimulate brain circuits and modulate neural activity.^{23} In our model, ultrasound frequency (*f* = 0.2 MHz) also adheres to low-frequency ultrasound standards. Overall, our TMAS strategy is capable of safely and reliably modulating epilepsy.

There are two potential limitations in our model. First, as observed in Fig. 5, ion concentrations such as $ [ K ] o$ and $ [ Na ] i$ do not return to the baseline level when neurons enter the quiescence state due to the termination of neuronal discharges or modulation by TMAS. We find out that this may be the consequence of the delayed action of the glial cells and excessive bath potassium concentration in our model. Another limitation is that we focus mainly on ultrasound’s temporal rather than spatial properties. Therefore, we did not consider the distribution of induced TMAS current in the brain. Our subsequent work may need to consider the effect of ultrasound attenuation on the skull and brain tissue.

Although our model focuses on a small-scale network of five neurons in the hippocampal CA3 region, we chose two classes of neurons including four pyramidal neurons and a PV+ interneuron, and complemented our neuron model with expressions for the ion concentration. First, the physiological ratio of pyramidal neurons to interneurons in the hippocampus is four to one. Second, some studies indicate that HYP seizure onset may depend on the preponderant involvement of principal (glutamatergic) neuron networks.^{7,8} Based on the two reasons mentioned above, we propose a quantitative configuration of these two types of neurons as a minimum model. In our paper, we concentrated on the transition of neuronal states induced by micro-level ionic concentration changes rather than the neuronal clustering activities in complex networks. Our results show HYP epilepsy evolution due to the depolarization block of PV+ interneurons and hyperexcitability of pyramidal neurons, which reproduce phenomena observed from physiological experiments. Inspired by *in vitro* experiments,^{50} our future work may construct a larger hippocampal CA3 region network where more GABAergic interneuron subtypes such as the cholecystokin interneuron and the somatostatin interneuron are added to study the effects of various interneurons on HYP seizures dynamics. In addition, another step may be to expand the model to include other brain regions to clarify the flow of information in pathological networks by nonlinear methods. Interestingly, according to Fig. 6, we find that the upper sine wave increases neuron firing rate while the lower sine wave decreases neuron firing rate. The different effects of ultrasound waveforms on neurons provide us with some inspiration. For neurological disorders such as epilepsy caused by neuronal over-excitation, we can choose the lower sine wave to inhibit neuron discharges. On the other hand, for neurological disorders such as major depressive disorder caused by neuronal over-inhibition, we can choose the upper sine wave to promote neuron discharges. Meanwhile, based on the high spatial resolution of TMAS, we can further explore potential target brain regions where TMAS can effectively modulate HYP seizure.

## ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (NNSFC) (Grant Nos. 12102014, 11932003, 32271361, and 12202022) and the National Key Research and Development Program of China (Grant Nos. 2021YFA1000200 and 2021YFC1000202). We thank the Figdraw platform (https://www.figdraw.com) for the visualization of Fig. 1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Zhiyuan Ma:** Conceptualization (lead); Data curation (lead); Writing – original draft (lead); Writing – review & editing (lead). **Yuejuan Xu:** Data curation (equal); Writing – review & editing (equal). **Gerold Baier:** Conceptualization (equal); Writing – review & editing (equal). **Youjun Liu:** Writing – review & editing (equal). **Bao Li:** Conceptualization (equal). **Liyuan Zhang:** Conceptualization (equal); Data curation (equal); Writing – original draft (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*in vitro*seizure-like events

^{2+}concentration with Chay neuron model

*Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting*

*in vitro*