Cortical spreading depression and spreading depolarization (CSD) are waves of neuronal depolarization that spread across the cortex, leading to a temporary saturation of brain activity. They are associated with various brain disorders such as migraine and ischemia. We consider a reduced version of a biophysical model of a neuron–astrocyte network for the initiation and propagation of CSD waves [Huguet *et al.,* Biophys. J. **111**(2), 452–462, 2016], consisting of reaction-diffusion equations. The reduced model considers only the dynamics of the neuronal and astrocytic membrane potentials and the extracellular potassium concentration, capturing the instigation process implicated in such waves. We present a computational and mathematical framework based on the parameterization method and singular perturbation theory to provide semi-analytical results on the existence of a wave solution and to compute it jointly with its velocity of propagation. The traveling wave solution can be seen as a heteroclinic connection of an associated system of ordinary differential equations with a slow–fast dynamics. The presence of distinct time scales within the system introduces numerical instabilities, which we successfully address through the identification of significant invariant manifolds and the implementation of the parameterization method. Our results provide a methodology that allows to identify efficiently and accurately the mechanisms responsible for the initiation of these waves and the wave propagation velocity.

We present a mathematical and computational analysis of the traveling wave solutions in a neuron–astrocyte model designed to capture the phenomenon of cortical spreading depolarization (CSD); brain waves of neuronal depolarization that propagate along the cortex. Our study introduces a powerful, highly efficient, and accurate methodology that yields reliable results for the computation of the traveling wave and its propagation velocity as well as a geometric understanding of the dynamics of CSD waves. We demonstrate that the method, based on the identification of significant invariant manifolds and the implementation of the parameterization method, can successfully overcome the numerical instabilities that arise from the inherent time scale differences in the system. Furthermore, our method paves the way for future research, enabling further investigations into the initiation and propagation mechanisms of CSD waves and facilitating exploration of more complex CSD wave models.

## I. INTRODUCTION

Cortical spreading depression and depolarization (CSD) are slowly propagating waves of rapid, near-complete depolarization of neurons and glial cells across the cortex.^{1–3} These waves move slowly through the cortex (2–6 mm/min), leading to a temporary cessation of brain electrical activity. This sustained depolarization of neurons has been implicated in various neurological disorders such as migraine and ischemia.^{2,4}

The hallmark of CSD is the disruption of transmembrane ionic gradients, which results in the collapse of ion homeostasis, swelling of cells, and elevated extracellular potassium levels.^{1,2,5} Astrocytes may play a crucial role in regulating potassium concentration by absorbing its accumulation from the extracellular space,^{6} thus potentially delaying or even preventing the onset of CSD waves.^{7–9} The ability of astrocytes to regulate potassium levels highlights the importance of astrocyte–neuron interaction in the control of brain function.

Dynamics involved in these phenomena occurs at different time scales. The processes governing ionic concentration changes occur at a slower time scale (seconds) compared to those that activate ionic currents and alter cell membrane potential (milliseconds).

In this work, we examine a biophysical model for cortical depolarization developed in Ref. 10. The model consists of a network of neurons and astrocytes, which was one of the first to incorporate a detailed description of astrocyte functioning. The interaction between neurons and astrocytes occurs through the diffusion of sodium and potassium ions in the extracellular space, as well as the release and uptake of these ions by the cells. This model was used in Ref. 10 to investigate, by means of numerical simulations, how the specific properties of astrocyte processes and the coupling between neurons and astrocytes, interact to generate or prevent the propagation of depolarization waves.

In this paper, we provide a reduction of the model in Ref. 10 based on the methodology developed in Ref. 11. The model is simplified so that it includes only those processes needed to initiate the wave. Thus, the reduced model considers only the dynamics of the neuronal and astrocytic membrane potentials and the extracellular potassium, while retaining the neuronal and glial pump currents, potassium release during firing, and diffusion along the extracellular space. This simplified framework makes the model more amenable to theoretical analysis.

We perform a study of the reduced model using computational and analytical tools to provide semi-analytical results on the existence of a traveling wave solution and its propagation velocity. Our approach consists in looking for a traveling wave solution by searching for a heteroclinic orbit in the resulting system of ordinary differential equations (ODEs). The presence of different times scales in the system allows us to apply slow–fast theory to study the system. Thus, we adopt the singular perturbation setting^{12,13} and introduce a small parameter that represents the ratio between the different time scales. Using techniques from classical singular perturbation theory,^{12,13} we break down the system into its fast and slow subsystems and construct a heteroclinic orbit in the singular limit. By leveraging the geometric understanding of the dynamics in the singular limit, we design a numerical strategy to compute the heteroclinic orbit beyond this limit, which effectively overcomes the numerical instabilities arising from the inherent time scale differences in the system. In our study, we present two distinct numerical strategies to achieve this, one based on the parameterization method^{14,15} and another employing Fenichel’s theory, which ensures the persistence of normally hyperbolic invariant manifolds.^{16–18}

Both methods have proven successful in mitigating the numerical instabilities associated with the system’s time scale differences. They yield accurate results for computing the traveling wave and its propagation velocity. Importantly, our methods shed light on the underlying mechanisms responsible for initiating the wave, providing valuable insights that can potentially be extended to the full model in the future.

The paper is organized as follows. In Sec. II, we present the original model introduced in Ref. 10 to describe CSD and the reduced version capturing the initiation of the wave. In Sec. III, we present the mathematical analysis of the traveling wave solutions of the reduced model using the singular perturbation framework and we construct a singular heteroclinic orbit. In Sec. IV, we provide the computation of the heteroclinic orbit and the velocity for the reduced system using the parameterization method.^{14,15} In Sec. V, we provide an alternative method to compute the heteroclinic orbit numerically based on Fenichel’s theory.^{16,17} We end with a discussion in Sec. VI. Finally, the Appendix contains some technical details involving algebraic mathematical computations.

## II. THE MATHEMATICAL MODEL

### A. Model for cortical spreading depression involving astrocytes

In this section, we present the neuron-astrocyte network model introduced in Ref. 10, upon which we base our work. The model describes quantitatively the dynamics of a network consisting of one array of neurons and another of astrocytes, sharing only the extracellular medium, through which there is diffusion of sodium (Na $ +$) and potassium (K $ +$) ions. The model captures accurately the initiation and propagation of cortical spreading depolarization (CSD) waves, emerging in several cortical diseases, such as migraine and ischemia.^{1} A thorough study on preventing or delaying this phenomenon was carried out in Ref. 10, with particular emphasis on the role played by the astrocyte cells and their gap junctions.

Parameters . | Paper . | Units . | Description . |
---|---|---|---|

R | 8.31 | J/mol K | Ideal gas constant |

F | 96 485 | C/mol | Faraday’s constant |

T | 310 | K | Absolute temperature |

C_{m} | 1 | μF/cm^{2} | Neuron’s membrane capacitance per unit area |

$ C m A$ | 1 | μF/cm^{2} | Astrocyte’s membrane capacitance per unit area |

ϕ_{n} | 0.8 | ms^{−1} | Maximum rate of activation of K^{+} channels |

$ \varphi h p$ | 0.05 | ms^{−1} | Maximum rate of activation of Na^{+} channels |

g_{Na} | 50.0 | mS/cm^{2} | Maximum fast Na^{+} current conductance |

g_{NaP} | 0.8 | mS/cm^{2} | Maximum fast Na^{+} current conductance |

g_{K} | 15 | mS/cm^{2} | Maximum delayed rectifier K^{+} conductance |

g_{L} | 0.5 | mS/cm^{2} | Nonspecific leak current conductance |

E_{L} | −70 | mV | Leak current reversal potential |

S_{N} | 922 | μm^{2} | Neuron membrane total surface area |

Ω_{n} | 2160 | μm^{3} | Neuron total volume |

S_{A} | 1600 | μm^{2} | Astrocyte membrane total surface area |

Ω_{a} | 2000 | μm^{3} | Astrocyte total volume |

Ω_{e} | α_{0} (Ω_{n} + Ω_{a}) | μm^{3} | Extracellular volume |

α_{0} | 0.2 | Volume fraction | |

D_{K} | 1.96 × 10^{−5} | cm^{2}/s | K^{+} Diffusion coefficient |

Δx | 0.044 | mm | Distance between neurons |

D_{Na} | 1.33 × 10^{−5} | cm^{2}/s | Na^{+} Diffusion coefficient |

P_{K} | 1 ⋅ 10^{−6} | cm/s | Membrane K^{+} permeability coefficient |

P_{Na} | 1.5 ⋅ 10^{−8} | cm/s | Membrane Na^{+} permeability coefficient |

ρ_{N} | 5 | μA/cm^{2} | Maximum neuron Na^{+}–K^{+} ATPase pump current |

ρ_{A} | 5 | μA/cm^{2} | Maximum astrocyte Na^{+}–K^{+} ATPase pump current |

Parameters . | Paper . | Units . | Description . |
---|---|---|---|

R | 8.31 | J/mol K | Ideal gas constant |

F | 96 485 | C/mol | Faraday’s constant |

T | 310 | K | Absolute temperature |

C_{m} | 1 | μF/cm^{2} | Neuron’s membrane capacitance per unit area |

$ C m A$ | 1 | μF/cm^{2} | Astrocyte’s membrane capacitance per unit area |

ϕ_{n} | 0.8 | ms^{−1} | Maximum rate of activation of K^{+} channels |

$ \varphi h p$ | 0.05 | ms^{−1} | Maximum rate of activation of Na^{+} channels |

g_{Na} | 50.0 | mS/cm^{2} | Maximum fast Na^{+} current conductance |

g_{NaP} | 0.8 | mS/cm^{2} | Maximum fast Na^{+} current conductance |

g_{K} | 15 | mS/cm^{2} | Maximum delayed rectifier K^{+} conductance |

g_{L} | 0.5 | mS/cm^{2} | Nonspecific leak current conductance |

E_{L} | −70 | mV | Leak current reversal potential |

S_{N} | 922 | μm^{2} | Neuron membrane total surface area |

Ω_{n} | 2160 | μm^{3} | Neuron total volume |

S_{A} | 1600 | μm^{2} | Astrocyte membrane total surface area |

Ω_{a} | 2000 | μm^{3} | Astrocyte total volume |

Ω_{e} | α_{0} (Ω_{n} + Ω_{a}) | μm^{3} | Extracellular volume |

α_{0} | 0.2 | Volume fraction | |

D_{K} | 1.96 × 10^{−5} | cm^{2}/s | K^{+} Diffusion coefficient |

Δx | 0.044 | mm | Distance between neurons |

D_{Na} | 1.33 × 10^{−5} | cm^{2}/s | Na^{+} Diffusion coefficient |

P_{K} | 1 ⋅ 10^{−6} | cm/s | Membrane K^{+} permeability coefficient |

P_{Na} | 1.5 ⋅ 10^{−8} | cm/s | Membrane Na^{+} permeability coefficient |

ρ_{N} | 5 | μA/cm^{2} | Maximum neuron Na^{+}–K^{+} ATPase pump current |

ρ_{A} | 5 | μA/cm^{2} | Maximum astrocyte Na^{+}–K^{+} ATPase pump current |

The currents $ I N a , gap$ and $ I K , gap$ in the astrocytes correspond to the potassium and sodium currents exchanged with nearby gap-junction-coupled astrocytes. They are modeled using the GHK formalism in Ref. 10. Since in this paper we will not consider them, we do not include details of the modelization here and we refer the reader to Ref. 10 for further information.

^{10,11}that increasing the pump strength, $ \rho N$ and $ \rho A$, increases the $ [ K + ] e$ threshold for wave initiation, which yields either a delay in the wave initiation (it takes longer to achieve the threshold) or a propagation failure (when the threshold is not achieved).

We performed a simulation of system (1) with the parameters given in Table I for a network of 50 neuron–astrocyte pairs sharing the extracellular space. To do so, we discretize the space variable $x$ using a grid of points $ x i=i\Delta x$ for $i=1,\u2026,N$, where $\Delta x$ is the distance between neurons. By employing finite differences and imposing Dirichlet boundary conditions at the end points of the array (which aims to simulate healthy tissue), we obtain a system of ODEs. See Appendix A for more details. We start the network at rest and inject an insult of $ K +$ to the extracellular space of the middle cells until a depolarizing wave is evoked. More precisely, the insult is modeled by adding a positive constant, $ I K , inj=0.005$ to the right-hand side of the equation for $ [ K + ] e$ corresponding to the middle cells until their voltage reaches the value $\u221230$ mV. The simulation shows a depolarization wave which propagates throughout the array of neuron–astrocyte pairs (see Fig. 1).

This solution corresponds to a traveling wave, a special type of solution of a partial differential equation (PDE) on an infinite domain with a strong restraint between $x$ and $t$ (i.e., $x\xb1ct$), whose spatial profile advances in time at a constant speed $c$ (therefore being invariant by temporal translations).

### B. Model reduction

In this section, we present a reduction of model (1) following the strategy employed in Ref. 11. First, we provide details of the simplifications and its main differences with respect to Ref. 11. Then, we present simulations of the reduced model, showing that it retains the ability to initiate and propagate a depolarizing front along the array of neuron–astrocyte pairs.

The strategy for the model reduction is to retain only the dynamics for the variables that play a significant role in the initiation of CSD waves, while neglecting the dynamics of the rest of variables. Thus, since the variation of sodium concentration is not essential for the wave initiation and propagation, the variables $ [ N a + ] e$, $ [ N a + ] i A$, and $ [ N a + ] i$ are frozen at the resting values. Namely, the resting values are set at 3.5 (mM) for the intracellular concentrations (i.e., $ [ N a + ] i= [ N a + ] i A=3.5$) and at 135 (mM) for the extracellular one (i.e., $ [ N a + ] e=135$).

The model is further simplified by removing the dynamics of the two gating variables $n$ and $ h p$. The variable $n$ has a small time constant, so it approaches fast the steady state $ n \u221e$, thus $n$ is set to $ n \u221e$. On the other hand, $ h p$ is a very slow variable, so we freeze it at its resting value. As we shall see later, the dynamics of $ h p$ plays a role in restoring the polarization of the neuron and astrocyte cells. When $ h p$ is held constant, the reduced system will show depolarizing fronts instead of pulses, which is sufficient to study the initiation and propagation of the waves.

We know from Refs. 10 and 7 that the astrocytes form a syncytium through gap junctions, which have a critical role in preventing wave propagation. Indeed, when the strength of the gap junctions is weak, a wave is initiated, as depicted in Fig. 1. However, in this paper, our focus is not on incorporating the effects of gap junctions. Instead, we intentionally remove the current through the astrocyte gap junctions $ I gap$ from the original model. This allows us to explore mathematical algorithms that capture the properties of the initiated wave.

In Fig. 2, we show a numerical simulation of the reduced model (8). As mentioned before, the resulting system can capture the initiation of the depolarization wave but fails to capture the cell’s membrane potential recovery phase to the resting potential, as the variable $ h p$ is frozen. Thus, it does no longer admits traveling pulse solutions (as displayed in Fig. 1) but traveling fronts. As in simulations of Fig. 1, we modeled an insult of $ K +$ in the extracellular space around the four middle neurons, in order to initiate the front.

From the numerical simulations, we can provide an estimate of the wave speed of the traveling front of system (8) by computing the elapsed time between the depolarization of two distinct cells. The speed is then obtained as the distance between the two cells divided by this time (see Fig. 3). We want to emphasize that we can rely on the numerical simulations to establish the existence of a traveling wave solution, but accurate estimates of its speed cannot be made due to the influence of the boundary conditions. Although we estimate the velocity using the voltage of cells that are far from the boundaries, we can still observe that the chosen boundary conditions (modeling healthy tissue at the boundaries of the domain) impact on the propagation speed of the traveling front by slowing it down. Indeed, our numerical explorations show that enlarging the size of the medium of the simulation by considering more neuron–astrocyte pairs, increases significantly the average wave speed (see Table II).

No. of neuron–astrocyte pairs . | Wave speed estimation (mm/min) . |
---|---|

50 | 3.8205 |

100 | 4.1641 |

300 | 4.7186 |

500 | 4.8253 |

No. of neuron–astrocyte pairs . | Wave speed estimation (mm/min) . |
---|---|

50 | 3.8205 |

100 | 4.1641 |

300 | 4.7186 |

500 | 4.8253 |

## III. TRAVELING WAVE SOLUTIONS IN THE SINGULAR LIMIT

. | $ p l 1$ . | $ p l 2$ . | p_{r}
. |
---|---|---|---|

x | −67.353 771 012 452 825 | −57.045 796 241 401 931 | 35.198 894 535 488 229 |

y | −63.416 145 863 486 385 | −55.561 014 704 831 557 | 11.631 018 842 324 311 |

z | 10.966 529 992 012 319 | 15.351 285 610 517 010 | 208.701 464 290 338 6 |

. | $ p l 1$ . | $ p l 2$ . | p_{r}
. |
---|---|---|---|

x | −67.353 771 012 452 825 | −57.045 796 241 401 931 | 35.198 894 535 488 229 |

y | −63.416 145 863 486 385 | −55.561 014 704 831 557 | 11.631 018 842 324 311 |

z | 10.966 529 992 012 319 | 15.351 285 610 517 010 | 208.701 464 290 338 6 |

By looking at the range of values taken by $ g ~$ in Fig. 4, one can argue that there exists an intermediate time scale between those of $x$ and $z$. This would give rise to another fast-slow version of system (13) with two parameters $\epsilon $ and $\delta $, which would have a term $\epsilon \delta $ multiplying $ x \u02d9$ (to indicate that the dynamics of $x$ is much faster than that of $y$). In this paper, though, we focus on the existence of two time scales and we leave this analysis for future work.

### A. Dynamics of the slow subsystem

*slow subsystem*—a system of differential-algebraic equations—given by

*critical manifold*$ M 0$, a two-dimensional manifold defined as

In Fig. 5(b), we show the manifold $ M ~ 0$ defined in (19), which is the projection of the critical manifold $ M 0$ defined in (18) onto the $(x,y,z)$ space (dashed gray curve) and a solution of the simulation of system (8) portrayed in Fig. 2 for a fixed value of the spatial variable (cyan curve, starting at the black dot). This illustrates that the solution stays close to the critical manifold $ M 0$ whenever possible, with a rapid transition between the lower branch $ M 0 l$ and the upper branch $ M 0 r$, that occurs at the fold line $ F R$. So, it is reasonable to start the study by computing an approximation of the heteroclinic orbit in the singular limit that consists of pieces of solutions of the slow and fast subsystems. This will be done in Sec. III C.

Notice that the equilibrium points for the slow subsystem coincide with the equilibrium points of the full system (13) even for $\epsilon \u22600$. The equilibrium points are given in Table III. Two of the equilibrium points, $ p l 1=( X l( z 1 l),Y( z 1 l), z 1 l,0)$ and $ p l 2=( X l( z 2 l),Y( z 2 l), z 2 l,0)$, lie on the lower branch $ M 0 l$ of the manifold and the third one, $ p r=( X r( z r),Y( z r), z r,0)$, on the upper branch $ M 0 r$.

### B. Dynamics of the fast subsystem

*fast subsystem*,

We will use the fast subsystem to study the stability of the critical manifold $ M 0$ (18) in its normal directions. Observe that, viewed in the fast subsystem, the critical manifold is actually a manifold filled with equilibrium points. The eigenvalues for the equilibrium points $( X l , m , r(z),Y(z),z,w)\u2208 R 4$ associated with the normal directions of the critical manifold are $ \lambda 1 \u2217= 1 c f x( X \u2217(z),z)$ with $\u2217=l,m,r$ and $ \lambda 2= 1 c g y(Y(z),z)<0$. In Fig. 8, we plot these eigenvalues as functions of $z$.

Observe that $ \lambda 2<0$ for all values of $z$. For points on $ M 0 l$, the eigenvalue $ \lambda 1 l<0$ and therefore, both eigenvalues are negative. The same happens for points on $ M 0 r$. However, for points on $ M 0 m$, $ \lambda 1 m>0$. From this analysis, we conclude that the lower and upper branches $ M 0 l$ and $ M 0 r$ are normally hyperbolic attracting manifolds, while the middle branch $ M 0 m$ is a normally hyperbolic saddle-type one.

On a final note, the stability of the fast subsystem’s equilibrium points (and hence that of the critical manifold’s branches) is reversed when the wave speed $c$ becomes negative: the lower and upper branches are then unstable and the middle one is of saddle type.

### C. Singular heteroclinic connection concatenating slow and fast subsystems

In Fig. 9 we plot, for different values of $c$, the unstable and stable invariant manifolds of $ p ^ l 1$ and $ p ^ r$, respectively (blue and orange curves). The invariant manifolds have been computed using the first-order approximation; we take as initial conditions $ p ^ \u2217+s v \u2217$ for $\u2217= l 1,r$, with $ v \u2217$ the eigenvectors of the linearization of the vector field (24) at $ p ^ \u2217$ and $s$ a small parameter, and integrate forward or backward in time accordingly. The (piecewise) heteroclinic orbit, plotted as solid black curve, is the solution of solving $d(c):= w u\u2212 w s=0$, where $ w u$ and $ w s$ denote the $w$-coordinates of the intersection points of the unstable and stable manifolds with the section $\Sigma $, respectively. The singular heteroclinic orbit is found for $ c 0\u2248 0.074 26$.

## IV. COMPUTING THE HETEROCLINC ORBIT FOR SYSTEM (15) WITH $\epsilon \u22600$ BY MEANS OF THE PARAMETERIZATION METHOD

In this section, we compute numerically the heteroclinic orbit of the full system (15) for $\epsilon \u22600$. Recall that for $\epsilon =0$ we had the critical normally hyperbolic attracting manifolds $ M 0 r$ and $ M 0 l$. Fenichel’s theory^{16,17} states that these critical manifolds persist, under small perturbations of $\epsilon $, to the full system (15) as $ M \epsilon l$ and $ M \epsilon r$, while preserving in turn their hyperbolicity properties. Moreover, the new invariant manifolds $ M \epsilon l$ and $ M \epsilon r$ contain the critical points $ p l 1$ and $ p r$, respectively.

Recall that manifolds $ M \epsilon l$ and $ M \epsilon r$ attract strongly all nearby orbits because of their normal directions being fast and contractive. Therefore, both equilibrium points $ p l 1$ and $ p r$ maintain their saddle-type nature but with one expanding and three contractive directions. This ensures that the heteroclinic orbit lives within the 1D unstable manifold of the point $ p l 1$, named $ \Gamma u$, until it arrives to the right fold, where it jumps to the upper branch and asymptotically approaches the manifold $ M \epsilon r$ and follows the 1D submanifold $ \Gamma s\u2282 W s( p r)$, tangent to the slowest contracting direction at $ p r$ (which is exponentially close to the unique stable direction of $ p r$ when restricted to the critical manifold $ M \epsilon r$).

Typically, detecting numerically a heteroclinic (or homoclinic) orbit involves following the unstable manifold of one equilibrium and the stable manifold of the other and determine whether they intersect on a convenient Poincaré section. Following manifolds of any dimension is however far from being a simple task. Even though in our present study the stable manifold of $ p r$ is three-dimensional, there is no need to follow the whole invariant object because we know the heteroclinic orbit will be through the 1D submanifold $ \Gamma s$.

Nevertheless, finding numerically a heteroclinic connection in the 4D system (15) proves to be a difficult task, due primarily to the presence of the two fast variables, $x$ and $y$, and the saddle character of the upper equilibrium point $ p r$. This prevents the usage of the same strategy used for the slow subsystem. On one side, following the 1D unstable manifold of $ p l 1$ $ \Gamma u$ works perfectly well until soon after the jump-up to the upper branch $ M \epsilon r$, where the orbit escapes along the unstable direction of the point $ p r$. On the other hand, integrating backwards $ \Gamma s$ blows up to infinity due to the presence of the two fast variables. Indeed, the fast variables compress the dynamics of the 4D space in the slow manifold very quickly, nonetheless backward integration reverses its stability, which turns $x$ and $y$ into expanding directions. Such fast expanding directions cause that orbits close to the 2D manifold $ M \epsilon r$ explode to infinity at the slightest numerical error.

For these reasons, following $ \Gamma s$ requires a more sophisticated method than simple backward integration from a linear approximation of the manifold. We will use the parameterization method^{15} in order to obtain a more accurate local approximation of $ \Gamma s$ which will be backward-integrated using an stiff solver (ode15s in Matlab).

The right-hand side of Eq. (34) can be numerically computed efficiently with accuracy as high as desired using the automatic differentiation algorithms.^{15,21} The algorithm outputs the coefficients (up to any order) of the power series of the composition of a truncated power series with any analytic function $F$. The method consists in computing systematically compositions of power series with elementary functions, for which there are well-established recurrences (see, for instance, Ref. 15).

Once fixed the order of the expansion $W(s)$ and given its coefficients $ W k$, we take the largest value $s$ such that the error in the invariance equation, $ | |F(W(s))\u2212DW(s)f(s) | | \u221e$, is less than $ 10 \u2212 10$ (see Fig. 10) and the backward integration from $W(s)$ provides a trajectory that hits the Poincaré section $ \Sigma 1:={z=22}$.

To obtain the 1D unstable manifold of the point $ p l 1$, $ \Gamma u$ we integrate the system forward starting from a point obtained using a first-order approximation of the manifold (i.e., $ p l 1+\delta v l 1$ with $ v l 1$ the eigenvector of the linearized system around $ p l 1$ associated to the positive eigenvalue and $\delta $ a small value).

In Figs. 11(a)–11(d), we illustrate the strategy used for computing a heteroclinic orbit. For a range of $c$-values including the value $ c 0= 0.074 26$ (which is the value for which we found a heteroclinic orbit in the singular case), we track both 1D invariant manifolds until they reach the Poincaré section $ \Sigma 1$ [see the 3D phase space projections in Figs. 11(a) and 11(c)]. The set of intersection points of the unstable manifold $ \Gamma u$ with the Poincaré section $ \Sigma 1$ forms a curve in terms of $c$ which intersects with that of the stable manifold $ \Gamma s$ [see Figs. 11(b) and 11(d)]. For a heteroclinic solution to exist, such an intersection should happen for the same value $c$. We therefore check that the distance $\Vert \u22c5 \Vert 2$ between the intersection points of both manifolds with the section $ \Sigma 1$ becomes zero for some value $ c ^$ [see the blue-dotted line in Fig. 11(f)]. In Fig. 11(f), the Euclidean distance has been disaggregated into the distances (with sign) of each of its components (see pink/purple/green-dotted lines). By applying a secant method in the $w$-component, we found an approximate value $ c ^=7.3135\xd7 10 \u2212 2$ with an error of order $ 10 \u2212 4$. This corresponds to an estimated wave velocity of $vel= c ^ D k=7.3135\xd7 10 \u2212 2 ms \u2212 1 / 2\u22c5 1.96 \xd7 10 \u2212 5 cm 2 s \u2212 1= 1.96\u22c5 0.073 135mm/s= 0.102 389mm/s=6.1433mm/min$.

## V. COMPUTING THE HETEROCLINIC ORBIT USING A DIFFERENT STRATEGY BASED ON FENICHEL’S THEORY

In this section, we present a different strategy to overcome the numerical difficulties to compute the 1D slow manifold $ \Gamma s$ of the critical point $ p r$. Our idea is to compute $ \Gamma s$ by integrating backwards system (15) starting from a local approximation of it, while restricting the trajectory to be approximately on $ M \epsilon r$, the perturbed upper branch of the manifold $ M 0 r$. By doing so, we assure that numerical integration will not blow up to infinity either in the $x$ or $y$ directions. As discussed in the previous section, the computation of the 1D unstable manifold of $ p l 1$, $ \Gamma u$, does not exhibit these problems and for this reason we focus only on the point $ p r$.

The coefficients will be determined by solving the invariance equation order by order. Actually, the coefficients of order 0 are given by $ m 0(z,w)= X r(z)$ and $ n 0(z,w)=Y(z)$ [see Eq. (20)].

The intersection points with section $ \Sigma 1$ will have $x$, $y$, and $w$ components for each value of $c$ [see Figs. 12(b) and 12(d)], and we look for $ c ~$ such that the triplet $ v u:=( x u, y u, w u)\u2208 \Gamma u$ equals $ v s:=( x s, y s, w s)\u2208 \Gamma s$. We compute the euclidean norm of the difference $ v u\u2212 v s$ and, separately, the distances with sign for each of its components as a function of the parameter $c$ [see Fig. 12(f)]. As we can observe, the euclidean distance vanishes within the range $(0.07,0.075)$, and using a secant method in the $w$-component we find that it vanishes for $ c ~= 0.073 135$ with an error of order $ 10 \u2212 4$. Notice that this value coincides with the one obtained using the parameterization method and corresponds to a wave velocity of $6.14mm/min$. In Fig. 12(e), we display the heteroclinic connection between $ p l 1$ and $ p r$ resulting for $ c ~$.

## VI. DISCUSSION

In this paper, we have presented a computational and mathematical framework to study the initiation and propagation of CSD waves in a reduction of the model introduced in Ref. 10, involving the dynamics of the extracellular potassium concentration $ [ K + ] e$. More precisely, we have studied semi-analytically the existence of traveling wave solutions in model (8), which correspond to heteroclinic orbits of the system (12). We have combined techniques from singular perturbation theory,^{12,13} and numerically efficient methods based on the parameterization method,^{15} involving automatic differentiation techniques, which provide an efficient computation of the heteroclinic solution and, subsequently, the propagation velocity of the wave. The methods have proven successful to overcome the difficulties associated with the presence of different time scales in the system.

We emphasize that the computation of the propagation velocity $vel\u2248c D K$ by means of computing a heteroclinic connection is more accurate since it does not have the boundary effects that appear in the simulation of the discretized PDE model (8) when the number of neurons is small. Indeed, the estimated velocity from the heteroclinic is faster (approx 6 mm/s) than the one obtained from numerical simulations of the PDE (approx 4 mm/s). Moreover, it is more computationally efficient since avoiding the boundary effects in the full model requires to use a large number of neurons, which, in turn, increases the computational cost of the simulations.

Our approach is based on the model reduction suggested in Ref. 11, However, we want to highlight some relevant differences with respect to this work. In Ref. 11, the system is further reduced to a single PDE for $ [ K + ] e$ of the reaction-diffusion type, imposing that the dynamics for $ V N$ and $ V A$ in (8) is set to be instantaneous (i.e., $d V n/dt=d V A/dt=0$). When imposing a traveling wave solution in this 1D reaction-diffusion PDE, it yields the dynamics described by the piecewise slow subsystem (24). In our study, we go beyond the singular limit and compute the heteroclinic orbit and the velocity in the case $\epsilon \u22600$, which corresponds to the full system (8). We emphasize though that the singular heteroclinic orbit provides a good approximation of the dynamics since the estimated velocity is similar and the trajectory tracks closely to the singular heteroclinic orbit [see black-dashed curves in Figs. 11(a) and 11(c) and 12(a) and 12(c)]. Indeed, in Fig. 13, we present the numerical simulations of (8) taking $ V N$ and $ V A$ as instantaneous variables (i.e., $ d V N d t= d V A d t=0$). As expected from our results, these simulations reveal good agreement with those for system (8) shown in Fig. 2.

The analysis performed herein was based on the existence of two different temporal scales in the reduced model (8): membrane potential variables change at a faster rate than the extracellular $ K +$ concentration. In the light of the observations in Fig. 4 (see Remark III.1), it seems reasonable to ponder another scale dynamics disparity, between the neuronal and astrocytic membrane potentials, and thus introduce a new intermediate scale $\delta $ on system (15). We leave this more detailed analysis for future work.

The presence of different time scales in a dynamical system can present a significant computational challenge. In this paper, we demonstrate the power and versatility of the parameterization method,^{14,15} widely used in the field of dynamical systems,^{22,23} to overcome some issues related to this challenge. Specifically, we applied the parameterization method to numerically compute the upper branch of the heteroclinic orbit in the system $\epsilon \u22600$ where conventional numerical methods failed due to the presence of a strong expansive directions.

We acknowledge that, even if our results on the existence of a traveling wave solution for the system with $\epsilon \u22600$ rely on accurate numerical computations, they are not a rigorous analytical proof of their existence. To do so, one needs to use techniques from singular perturbation theory^{12,13} to describe the behavior of the solutions close to the fold of the slow manifold and techniques from computer assisted proofs^{24} to keep track of the stable and unstable manifolds of the saddle points, until they reach suitable Poincaré sections, which we plan to attempt in future studies.

In this paper, we have assumed a reduced model (8) without gap junctions in the astrocyte cells. However, previous work^{10} has shown that gap junctions play a significant role in preventing wave propagation. Therefore, future studies will involve incorporating gap junction currents and exploring how they affect the wave shape and dynamics.

Finally, we want to emphasize that the reduced model (8) only includes the biophysics terms affecting the wave initiation but does not model the wave termination. Indeed, $N a +$ concentrations are assumed constant in the reduced model, but in the original model presented in Ref. 10, the increase in $ [ K + ] e$ is accompanied by a steep decline in $ [ N a + ] e$ (see Fig. 1). Moreover, freezing the inactivation variable $ h p$ of the persistent sodium current $ I N a P$ prevents the termination of the wave and retains the network cells in the depolarized state. Further analysis may focus on the impact of including the dynamics for these variables, by designing a different reduced system which, at this turn, involves only the biophysics processes involved in the wave termination. The description of the wave termination would be analogous to that conducted within this manuscript, but the heteroclinic orbit would go from the equilibrium point corresponding to the depolarized state to the one corresponding to the resting state. Thus, in the future, we plan to concatenate two reduced systems, one that describes the wave initiation and another describing the wave termination, using a piecewise formulation.^{20,25,26} This more accurate description of the wave will provide new insights into the dynamics of the system describing CSD, leading to a more comprehensive understanding of the mechanisms underlying the CSD phenomenon.

In summary, we have presented a computational and mathematical framework that offers a more computationally efficient and accurate approach to identify the properties of CSD waves in model (1) compared to the crude simulation of the model. We hope that it establishes a solid foundation upon which future studies of more complex models can be built (such as extensions of the models in Ref. 27 incorporating the spatial domain).

## ACKNOWLEDGMENTS

Work produced with support of the grant PID-2021-122954NB-100 funded by MCIN/AEI/ 10.13039/501100011033 and “ERDF: A way of making Europe.” T.M.S and G.H acknowledge the Maria de Maeztu Award for Centers and Units of Excellence in R&D (No. CEX2020-001084-M). T.M.S. is supported by the Catalan Institution for Research and Advanced Studies via an ICREA Academia Prize 2019. We also acknowledge the use of the UPC Dynamical Systems group’s cluster for research computing (https://dynamicalsystems.upc.edu/en/computing/)

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**David Reyner-Parra:** Conceptualization (lead); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Software (lead); Supervision (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). **Carles Bonet:** Conceptualization (supporting); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (supporting). **Teresa M. Seara:** Conceptualization (supporting); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Gemma Huguet:** Conceptualization (lead); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Software (lead); Supervision (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

### APPENDIXES

In this section we include details of the numerical and algebraic mathematical computations.

### APPENDIX A: SOLVING NUMERICALLY THE PDE SYSTEMS

In this section, we discuss the numerical details for the simulations of the PDE systems (1), (8), and (8) with the reductions in (9) (Figs. 1, 2, and 13, respectively). We use a network of N neuron–astrocyte pairs sharing the same extracellular space. Thus, the space is discretized taking a grid $ x i=i\Delta x$, for $i=1,\u2026,N$ and $\Delta x$ is the distance between neurons.

The initial conditions are identical for all the pairs. We used the stiff integrator ode15s of Matlab with crude tolerances.

### APPENDIX B: COMPUTATION OF THE DERIVATIVES OF THE FUNCTIONS $ H \u2217$

### APPENDIX C: TAYLOR EXPANSIONS OF THE MANIFOLD $ M \epsilon r$ in $\epsilon $

To compute Taylor expansions of the manifold $ M \epsilon r$ in $\epsilon $, we solve the invariance equation (36) by matching powers of $\epsilon $.

After substituting the Taylor expansions (C1) into system (36), we are able to solve the resulting system of equations up to any power of $\epsilon $,

- $O(1)$:$f(m0,z)=0\u27f9m0(z,w)=Xr(z)(see the red curve in {Fig. 5}),g(n0,z)=0\u27f9n0(z,w)=Y(z)(see the orange curve in {Fig. 5}).$
- $O(\epsilon )$:Observe that, since$f(Xr(z),z)=0$and $g(Y(z),z)=0$, we have used(C2)$\u2202m0\u2202z=dXr(z)dz=\u2212fz(m0,z)fx(m0,z)and\u2202n0\u2202z=dYdz=\u2212gz(n0,z)gy(n0,z).$
- $O(\epsilon 2)$ :

### APPENDIX D: COMPUTATION OF THE PARTIAL DERIVATIVES OF THE FUNCTIONS $ f ~(x,z)$ AND $ g ~(y,z)$

### APPENDIX E: LINEAR LOCAL APPROXIMATION FOR THE 1D STABLE MANIFOLD $ \Gamma ~ s$

## REFERENCES

*Dynamical Systems. Lecture Notes in Mathematics*, edited by R. Johnson (Springer, Berlin, 1995), Vol. 1609.

*Multiple Time Scale Dynamics. Applied Mathematical Sciences*(Springer, Cham, 2015), Vol. 191.

*The Parameterization Method for Invariant Manifolds*

*Lecture Notes in Mathematics*(Springer, Berlin, 1977), Vol. 538.

*Mathematical Foundations of Neuroscience*

*Interdisciplinary Applied Mathematics*(Springer, 2008).

*Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation*