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.

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 setting12,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 method14,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.

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.

The mathematical model consists of a set of ordinary and partial differential equations modeling the dynamics of the neuronal and astrocytic membrane potentials, V N and V A, respectively, the extracellular and intracellular concentrations of sodium and potassium, [ N a + ] e, [ N a + ] i, [ N a + ] i A, [ K + ] e, [ K + ] i, [ K + ] e A, respectively, and the gating variables n and h p. The system of differential equations is given by
C m d V N d t = I N a I N a P I K I L I P m , C m A d V A d t = I N a A I K A I P m A I gap , d n d t = ϕ n n ( V N ) n τ n ( V N ) , d h p d t = ϕ h p h p ( V N ) h p τ h p ( V N ) , d [ N a + ] i d t = 10 S N F Ω n ( I N a + I N a P + 3 I P m ) , d [ N a + ] i A d t = 10 S A F Ω a ( I N a A + 3 I P m A I N a , gap ) , d [ K + ] i d t = 10 S N F Ω n ( I K 2 I P m ) , d [ K + ] i A d t = 10 S A F Ω a ( I K A 2 I P m A I K , gap ) , [ N a + ] e t = D N a 2 [ N a + ] e x 2 + 10 S N F Ω e ( I N a + I N a P + 3 I P m ) + 10 S A F Ω e ( I N a A + 3 I P m A ) , [ K + ] e t = D K 2 [ K + ] e x 2 + 10 S N F Ω e ( I K 2 I P m ) + 10 S A F Ω e ( I K A 2 I P m A ) .
(1)
The ionic currents I N a, I N a P, I K, and I L are modeled as in the Hodgkin–Huxley model, that is,
I N a = g N a m 3 ( V N ) ( 1 n ) ( V N E N a ) , I N a P = g N a P m p 3 ( V N ) h p ( V N E N a ) , I K = g K n 4 ( V N E K ) , I L = g L ( V N E L ) .
(2)
The reversal potential E X, for X { N a , K }, is modeled by the Nernst equation
E X = R T F ln [ X ] e [ X ] i ,
(3)
where R , T, and F are the gas constant, temperature, and Faraday’s constant, respectively. The values are given in Table I.
The asymptotic functions m , m p , n , and h p modeling the steady state values of the ionic gating variables are sigmoidal-shaped functions of the form
X ( V ) = 1 1 + e ( ( V V X ) / θ X ) ,
(4)
for X { m , n , m p , h p }, with V m = 34, θ m = 5, V n = 55, θ n = 14, V m p = 40, θ m p = 6, V h p = 48, and θ h p = 6. Finally, the time constants modeling the rate at which the gating variables reach its steady state values are
τ n ( V ) = 0.05 + 0.27 1 + e ( ( V + 40 ) / ( 12 ) ) , τ h p ( V ) = 10000 cosh ( ( V + 48 ) / 12 ) .
(5)
The ionic currents I N a A and I K A in the modelization of the astrocyte equation are described by the Goldman–Hodgkin–Katz (GHK) formulation,
I N a A = P N a F ϕ [ N a + ] e e ϕ [ N a + ] i A e ϕ 1 , I K A = P K F ϕ [ K + ] e e ϕ [ K + ] i A e ϕ 1 ,
(6)
with ϕ := F R T V A, P N a = 1.5 × 10 8 and P K = 1 × 10 6. Notice that the currents I N a A and I K A in (6) are not defined when the astrocyte’s membrane potential reaches zero because the dividing term e ϕ 1 cancels. This is resolved by explicitly defining the values of I N a A and I K A at V A = 0,
I N a A ( V A , [ N a + ] e , [ N a + ] i A ) | V A = 0 := lim V A 0 P N a F ϕ [ N a + ] e e ϕ [ N a + ] i A e ϕ 1 = P N a F ( [ N a + ] i A [ N a + ] e ) , I K A ( V A , [ K + ] e , [ K + ] i A ) | V A = 0 := lim V A 0 P K F ϕ [ K + ] e e ϕ [ K + ] i A e ϕ 1 = P K F ( [ K + ] i A [ K + ] e ) .

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.

The sodium–potassium pump currents in neurons and astrocytes, I P m and I P m A, respectively, are given by
I P m = ρ N ( [ K + ] e 2 + [ K + ] e ) 2 ( [ N a + ] i 7.7 + [ N a + ] i ) 3 , I P m A = ρ A ( [ K + ] e 2 + [ K + ] e ) 2 ( [ N a + ] i A 7.7 + [ N a + ] i A ) 3 ,
(7)
where ρ N and ρ A represent the maximal current allowed to go through the neuron and astrocyte pumps, respectively (regarded as parameters in Ref. 10 to explore different conditions). In the present work, we have fixed both of them at a constant value of 5, for which we know from Ref. 10 that a wave can be initiated and propagated through the tissue. We know from previous works10,11 that increasing the pump strength, ρ N and ρ 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 Δ x for i = 1 , , N, where Δ 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 30 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 ± c t), whose spatial profile advances in time at a constant speed c (therefore being invariant by temporal translations).

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).

On the other hand, since the astrocytes’ volume is bigger than that of the extracellular medium [i.e., Ω e = α 0 ( Ω n + Ω a ), α 0 1], changes in the quantity of K + ions will have a greater impact on the extracellular K + concentration [ K + ] e than on the concentration within the astrocytes [ K + ] i A. Thus, the variable [ K + ] i A in Eq. (1) is also taken as constant. Finally, the conservation of K + ions states that
[ K + ] i := 1 Ω n ( K tot + Ω e [ K + ] e Ω a [ K + ] i A ) ,
being K tot + the total initial number of K + ions in the space. Note that, unlike [ N a + ] i, the magnitude [ K + ] i is not constant as it depends moderately on the slowly changing dynamics of [ K + ] e. The initial values for the extracellular and intracellular K + concentrations are [ K + ] e = 3.5 (mM) and [ K + ] i = 135 (mM), respectively, while [ K + ] i A is kept constant at 135 (mM).

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 , thus n is set to n . 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.

These simplifications render system (1) depending only on the variables V N, V A, and [ K + ] e in the following way:
C m d V N d t = I N a I N a P I K I L I P m , C m A d V A d t = I N a A I K A I P m A , [ K + ] e t = D K 2 [ K + ] e x 2 + 10 S N F Ω e ( I K 2 I P m ) + 10 S A F Ω e ( I K A 2 I P m A ) .
(8)
System (8) is the model we will analyze in this paper, which is the same considered in Ref. 11 except that we keep the terms I K A and I P m A in the equation for [ K + ] e, describing the influence of astrocytes in the extracellular levels of K +.

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).

A further simplification of model (8) has been done in Ref. 11, by leveraging the different scales in the rate of change of the membrane potential variables, V N and V A, and the extracellular K + (compare, for instance, the evolution of V N and V A relative to [ K + ] e in Fig. 2). Mathematically, imposing that
I N a I N a P I K I L I P m = 0 , I N a A I K A I P m A = 0 ,
(9)
the system is reduced to a single PDE for [ K + ] e of the reaction-diffusion type, that preserves qualitatively the dynamics shown in Fig. 2. As we will see later, part of our study of system (8) will require the knowledge of the dynamics of this further reduced model.
In this section, we present a strategy that combines analytical and numerical methods to show that system (8) has traveling wave solutions and to compute them. The first step is to postulate a solution of the form
V N ( x , t ) = V ¯ N ( ξ ) , V A ( x , t ) = V ¯ A ( ξ ) , with ξ := x / D K + c t , [ K + ] e ( x , t ) = [ K + ] ¯ e ( ξ ) ,
(10)
where the constant c (to be determined) establishes the wave propagation velocity, v e l = c D K. Substituting (10) into (8) results in a system of ordinary differential equations, given by
c d V ¯ N d ξ = 1 C m ( I N a + I N a P + I K + I L + I P m ) , c d V ¯ A d ξ = 1 C m A ( I N a A + I K A + I P m A ) , c d [ K + ] ¯ e d ξ = d 2 [ K + ] ¯ e d ξ 2 + 10 S N F Ω e ( I K 2 I P m ) + 10 S A F Ω e ( I K A 2 I P m A ) .
(11)
We will see that system (11) has two saddle points p l 1 and p r (see Table III). To find a traveling front solution of system (8), we will look for a heteroclinic orbit of system (11) from p l 1 to p r. The goal is to determine whether the parameter c can be chosen so that this heteroclinic orbit exists. We will do the analysis for c > 0.
We relabel the variables V ¯ N, V ¯ A, and [ K + ] ¯ e as x (not to be confused with the spatial variable x), y and z, respectively, and introduce the variable w := d [ K + ] ¯ e d ξ = d z d ξ, to finally obtain
c d x d ξ = 1 C m ( I N a + I N a P + I K + I L + I P m ) , c d y d ξ = 1 C m A ( I N a A + I K A + I P m A ) , d z d ξ = w , d w d ξ = c w 10 S N F Ω e ( I K 2 I P m ) 10 S A F Ω e ( I K A 2 I P m A ) ,
(12)
which, for simplicity, can be written in a more general form,
d x d ξ = 1 c f ~ ( x , z ) , d y d ξ = 1 c g ~ ( y , z ) , d z d ξ = w , d w d ξ = c w h ( x , y , z ) ,
(13)
with
f ~ ( x , z ) = 1 C m ( I N a ( x ) + I N a P ( x ) + I K ( x , z ) + I L ( x ) + I P m ( z ) ) , g ~ ( y , z ) = 1 C m A ( I N a A ( y ) + I K A ( y , z ) + I P m A ( z ) ) , h ( x , y , z ) = 10 S N F Ω e ( I K ( x , z ) 2 I P m ( z ) ) + 10 S A F Ω e ( I K A ( y , z ) 2 I P m A ( z ) ) ,
(14)
where we have explicitly indicated the dependence on x, y, and z for each ionic current.
We observe that the neuronal and astrocytic membrane potentials ( x and y variables, respectively) change more rapidly than that of the extracellular concentration of potassium in the medium ( z variable). In Fig. 4, we plot functions f ~ ( x , z ), g ~ ( y , z ), and h ( x , y , z ) to illustrate so. Indeed, the ranges of f ~ and g ~ are larger than that of h, entailing a higher rate of change in the derivatives d x d ξ and d y d ξ and thus a much faster dynamics. For this reason, we consider system (13) with an artificial small parameter ε to emphasize these slow–fast dynamics, that is,
ε d x d ξ = 1 c f ( x , z ) , ε d y d ξ = 1 c g ( y , z ) , d z d ξ = w , d w d ξ = c w h ( x , y , z ) ,
(15)
where we have denoted f = ε f ~ and g = ε g ~. Equivalently, re-scaling the time variable ε η = ξ, system (15) becomes
d x d η = 1 c f ( x , z ) , d y d η = 1 c g ( y , z ) , d z d η = ε w , d w d η = ε ( c w h ( x , y , z ) ) .
(16)
Remark III.1

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 ε and δ, which would have a term ε δ multiplying x ˙ (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.

Note that Eqs. (15) and (16) correspond to a singular perturbation problem in parameter ε. First, we are going to construct a heteroclinic orbit in the singular limit, that is when ε 0 by splitting the dynamics into the slow and fast dynamics.

Taking ε = 0 in system (15), we obtain the so-called slow subsystem—a system of differential-algebraic equations—given by
0 = 1 c f ( x , z ) ,
(17a)
0 = 1 c g ( y , z ) ,
(17b)
d z d ξ = w ,
(17c)
d w d ξ = c w h ( x , y , z ) .
(17d)
The first two algebraic equations (17a) and (17b) describe the commonly named critical manifold M 0, a two-dimensional manifold defined as
M 0 := { ( x , y , z , w ) R 4 | f ( x , z ) = 0 , g ( y , z ) = 0 } := M ~ 0 × R ,
(18)
where
M ~ 0 := { ( x , y , z ) R 3 | f ( x , z ) = 0 , g ( y , z ) = 0 } .
(19)
Figure 5(a) shows the curves f ( x , z ) = 0 and g ( y , z ) = 0 on the ( x , z ) and ( y , z )-plane, respectively. Notice that for a finite range of values of z, the equation f ( x , z ) = 0 has three branches x = X l , m , r ( z ), where when they are defined,
X l ( z ) X m ( z ) X r ( z ) .
If we denote by ( x L , z L ) and ( x R , z R ), the two fold points [i.e., solutions of the system f ( x , z ) = 0, f x ( x , z ) = 0], such that z L < z R, then X l is defined for z < z R, X m for z L < z < z R and X r for z > z L. Notice though that the curve g ( y , z ) = 0 satisfies g y ( y , z ) < 0 for all z and therefore defines y = Y ( z ) globally [see Fig. 5(a)]. Thus, the manifold M 0 has three branches [see Fig. 5(b)], namely,
M 0 l = { ( x , y , z , w ) R 4 | z < z R , x = X l ( z ) , y = Y ( z ) } , M 0 m = { ( x , y , z , w ) R 4 | z L < z < z R , x = X m ( z ) , y = Y ( z ) } , M 0 r = { ( x , y , z , w ) R 4 | z > z L , x = X r ( z ) , y = Y ( z ) } ,
(20)
as well as two fold lines F L , R = { ( x L , R , Y ( z L , R ) , z L , R ) } × R.

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.

The next step is to study the dynamics of the slow subsystem (17), which describes the motion on M 0. The dynamics on the branch M 0 for = l , m , r, of the manifold M 0 is given by system
{ d z d ξ = w , d w d ξ = c w h ( X ( z ) , Y ( z ) , z ) ,
(21)
where X ( z ) and Y ( z ) are the functions defining the manifolds M 0 introduced in (20).
We begin by looking at the equilibrium points of system (21). Therefore, we seek for zeros of
H ( z ) := h ( X ( z ) , Y ( z ) , z ) ,
(22)
for = l , m , r. Notice that each function H is defined for a different range of z. In Fig. 6, we plot the functions H ( z ) (distinguished by the color), each one defined in its respective domain. It turns out that H l has two zeros, z 1 l and z 2 l, H r only one, z r, while H m none.

Notice that the equilibrium points for the slow subsystem coincide with the equilibrium points of the full system (13) even for ε 0. 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.

If we denote by F s s ( z , w ), for = l , m , r, the vector fields of system (21), the stability of the equilibrium points p is determined by the eigenvalues λ ± of the matrix
D F s s ( z , 0 ) = ( 0 1 d H d z c ) .
Observe that det D F s s = d H d z and Tr D F s s = c > 0. The expressions for the derivatives of H are given in  Appendix B. Function H l has a negative slope at point z 1 l and so does H r at z r, hence they are saddle points. The derivative of H l at z 2 l is positive, and thereby assuring its unstable character. To illustrate this, in Fig. 7 we plot both eigenvalues λ ± for each equilibrium point for a range of positive values c.
Setting now ε = 0 in system (16), we obtain the so-called fast subsystem,
d x d η = 1 c f ( x , z ) ,
(23a)
d y d η = 1 c g ( y , z ) ,
(23b)
d z d η = 0 ,
(23c)
d w d η = 0 .
(23d)

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 ) R 4 associated with the normal directions of the critical manifold are λ 1 = 1 c f x ( X ( z ) , z ) with = l , m , r and λ 2 = 1 c g y ( Y ( z ) , z ) < 0. In Fig. 8, we plot these eigenvalues as functions of z.

Observe that λ 2 < 0 for all values of z. For points on M 0 l, the eigenvalue λ 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, λ 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.

In this section, we look for a heteroclinic connection between points p l 1 and p r in the singular limit combining solutions from the slow and fast subsytems (17) and (23). Since we are looking for an orbit that tracks the critical manifold branch M 0 l until it reaches the fold F R, where it rapidly switches to track the upper branch M 0 r, we build a piecewise system defined on the ( z , w ) plane of the form
d z d ξ = w , d w d ξ = c w H ( z ) ,
(24)
where H ( z ) is a piecewise smooth function defined as
H ( z ) = { H l ( z ) = h ( X l ( z ) , Y ( z ) , z ) if z z R , H r ( z ) = h ( X r ( z ) , Y ( z ) , z ) if z z R ,
where z R = 18.276 is the z-coordinate of the right fold F R.
We look for a heteroclinic orbit between the two saddle points, p ^ l 1 = ( z 1 l , 0 ) and p ^ r = ( z r , 0 ), in the two-dimensional (piecewise) system (24). To that end, we must compute the unstable and stable invariant manifolds of p ^ l 1 and p ^ r, respectively, which are 1D curves, and extend them until reaching the section
Σ := { ( z , w ) R 2 | z = z R } ,
(25)
where the definition of function H ( z ) changes from the lower branch M 0 l to the upper branch M 0 r of the critical manifold. Generically, these pieces of curves will not meet Σ at the same point. For this to occur, we treat c as a parameter and seek for a value c 0 such that the distance between the intersection points of both invariant manifolds on Σ is zero. Notice that the resulting heteroclinic orbit will be piecewise-defined.

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 ^ + s v for = l 1 , r, with v the eigenvectors of the linearization of the vector field (24) at p ^ 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 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 Σ, respectively. The singular heteroclinic orbit is found for c 0 0.074 26.

We recognize (24) as the bistable reaction-diffusion equation studied in Ref. 11. Adapting classical methods described in Refs. 20 and 19, one can show that the function H ( z ) satisfies the hypothesis that guarantees the existence of heteroclinic orbits for system (24).

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

Recall that manifolds M ε l and M ε 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 Γ u, until it arrives to the right fold, where it jumps to the upper branch and asymptotically approaches the manifold M ε r and follows the 1D submanifold Γ s 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 ε 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 Γ 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 Γ u works perfectly well until soon after the jump-up to the upper branch M ε r, where the orbit escapes along the unstable direction of the point p r. On the other hand, integrating backwards Γ 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 ε r explode to infinity at the slightest numerical error.

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

Next, we provide details of the parameterization method applied to our problem. We seek for a parametrization of the manifold Γ s as
W : R R 4 s W ( s ) .
(26)
Simultaneously, we will also determine the internal dynamics of the manifold, governed by the ODE s ˙ = f ( s ) with f : R R. Imposing W ( s ) to be a solution of the system X ˙ = F ( X ) (system (15) in our case), we obtain the invariance equation
F ( W ( s ) ) = D W ( s ) f ( s ) .
(27)
Assume W ( s ) and f ( s ) can be both expressed as power series in s,
W ( s ) = k = 0 W k s k , f ( s ) = k = 0 f k s k ,
(28)
where coefficients W k R 4 and f k R. The first coefficient of the parametrization, W 0, is actually the equilibrium point (in our case p r). After substituting (28), the invariance equation (27) is left as
F ( k = 0 W k s k ) = ( k = 1 W k k s k 1 ) ( k = 0 f k s k ) ,
to be solved order-by-order to determine the unknown coefficients W k and f k. One easily obtains that f 0 = 0 so that equation of order 0 holds. For higher orders, we develop the left-hand side of the above equation by Taylor expanding the vector field F about the equilibrium point
F ( p r ) + D F ( p r ) ( k 1 W k s k ) + 1 2 ! D 2 F ( p r ) ( k 1 W k s k ) 2 + k 0 F k s k = ( k = 1 W k k s k 1 ) ( k = 1 f k s k ) .
(29)
The equation of order 1 leads to the eigenvalue problem
D F ( p r ) W 1 = f 1 W 1 ,
(30)
thus being W 1 the eigenvector of eigenvalue f 1. Since we are looking for the slow (sub)manifold (associated to the slowest eigenvalue), we choose f 1 := λ slow and W 1 the corresponding eigenvector, that we will take of modulus 1.
Solving the equations of higher order shows that there is not a unique solution. We will use this freedom to choose the dynamics f to be as simple as possible, namely, as
f ( s ) = λ slow s ,
(31)
which corresponds to f k = 0 for k 2.
With this simplification, the equation of order 2 reads as
D F ( p r ) W 2 + [ F ( W < 2 ) ] | 2 F 2 = 2 λ slow W 2 ,
(32)
where W < 2 = k < 2 W k s k and [ F ( W < 2 ) ] | 2 denotes the coefficient of order 2 of the Taylor expansion in s. The resulting system is linear
( D F ( p r ) 2 λ slow Id ) W 2 = [ F ( W < 2 ) ] | 2 ,
(33)
and, therefore, it can be solved efficiently. In general, for the coefficient of order k 2, W k, we need to solve the linear system
( D F ( p r ) k λ slow Id ) W k = [ F ( W < k ) ] | k ,
(34)
where, again, the coefficient F k has been split into the known, [ F ( W < k ) ] | k, and unknown parts, D F ( p r ) W k.
Equation (34) can be solved provided
det ( D F ( p r ) k λ slow Id ) 0 ,
or, equivalently, that k λ slow is not an eigenvalue of the differential matrix D F ( p r ), which is fulfilled in our case because the other eigenvalues are much larger (this is equivalent to say there are no resonances among the eigenvalues of the system).

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 ) ) D W ( s ) f ( s ) | | , is less than 10 10 (see Fig. 10) and the backward integration from W ( s ) provides a trajectory that hits the Poincaré section Σ 1 := { z = 22 }.

To obtain the 1D unstable manifold of the point p l 1, Γ u we integrate the system forward starting from a point obtained using a first-order approximation of the manifold (i.e., p l 1 + δ v l 1 with v l 1 the eigenvector of the linearized system around p l 1 associated to the positive eigenvalue and δ 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 Σ 1 [see the 3D phase space projections in Figs. 11(a) and 11(c)]. The set of intersection points of the unstable manifold Γ u with the Poincaré section Σ 1 forms a curve in terms of c which intersects with that of the stable manifold Γ 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 2 between the intersection points of both manifolds with the section Σ 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 × 10 2 with an error of order 10 4. This corresponds to an estimated wave velocity of v e l = c ^ D k = 7.3135 × 10 2 ms 1 / 2 1.96 × 10 5 cm 2 s 1 = 1.96 0.073 135 mm/s = 0.102 389 mm/s = 6.1433 mm/min.

In this section, we present a different strategy to overcome the numerical difficulties to compute the 1D slow manifold Γ s of the critical point p r. Our idea is to compute Γ s by integrating backwards system (15) starting from a local approximation of it, while restricting the trajectory to be approximately on M ε 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, Γ u, does not exhibit these problems and for this reason we focus only on the point p r.

Since the upper branch of the Fenichel’s manifold M ε r perturbs from the critical manifold M 0 r, it can be approximated by considering both fast variables x and y as formal expansions in ε,
x = m ( z , w , ε ) = m 0 ( z , w ) + ε m 1 ( z , w ) + ε 2 m 2 ( z , w ) + , y = n ( z , w , ε ) = n 0 ( z , w ) + ε n 1 ( z , w ) + ε 2 n 2 ( z , w ) + .
(35)

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)].

From imposing the manifold (35) to be invariant by system (15), we reach the following system of equations:
x ˙ = m ˙ ( z , w , ε ) f ( m ( z , w , ε ) , z ) = ε c w m z + ε c ( c w h ( m ( z , w , ε ) , n ( z , w , ε ) , z ) ) m w , y ˙ = n ˙ ( z , w , ε ) g ( n ( z , w , ε ) , z ) = ε c w n z + ε c ( c w h ( m ( z , w , ε ) , n ( z , w , ε ) , z ) ) n w .
(36)
Developing functions f, g, and h as Taylor series around the coefficients of order 0, m 0 and n 0, and substituting them back into system (36), we are able to solve the resulting system of equations up to any power of ε. Indeed, one can obtain the following expressions (see  Appendix C) up to order 2:
m 1 ( z , w ) = c w f z ( X r ( z ) , z ) ( f x ( X r ( z ) , z ) ) 2 , n 1 ( z , w ) = c w g z ( Y ( z ) , z ) ( g y ( Y ( z ) , z ) ) 2 ,
and
m 2 ( z , w ) = 1 2 ! f x x ( m 0 , z ) f x ( m 0 , z ) m 1 2 + c w f x ( m 0 , z ) m 1 z c ( c w h ( m 0 , n 0 , z ) ) c f z ( m 0 , z ) ( f x ( m 0 , z ) ) 3 ,
and
n 2 ( z , w ) = 1 2 ! g y y ( n 0 , z ) g y ( n 0 , z ) n 1 2 + c w g y ( n 0 , z ) n 1 z c ( c w h ( m 0 , n 0 , z ) ) c g z ( n 0 , z ) ( g y ( n 0 , z ) ) 3 .
Next, we will restrict the dynamics of the full system (15) onto the second-order approximation of the manifold M ε r by taking
x = X r ( z ) + ε m 1 ( z , w ) + ε 2 m 2 ( z , w ) , y = Y ( z ) + ε n 1 ( z , w ) + ε 2 n 2 ( z , w ) .
(37)
The resulting system is
z ˙ = ε w , w ˙ = ε ( c w h ( X r ( z ) + ε m 1 ( z , w ) + ε 2 m 2 ( z , w ) , Y ( z ) + ε n 1 ( z , w ) + ε 2 n 2 ( z , w ) , z ) ) .
(38)
We compute the 1D stable manifold Γ ~ s of the critical point ( z r , 0 ) of system (38). We take the local approximation provided by the linear terms and we integrate it backwards in time until the trajectory hits the Poincaré section Σ 1 = { z = 22 }. Details are given in  Appendix E. Thus, Γ s can be obtained embedding Γ ~ s in R 4 using (37). In Figs. 12(a) and 12(c), we show manifolds Γ s and Γ u (where the latter has been computed as described in the previous section) and their intersection along the Poincaré section Σ 1 (depicted in dark gray), for different values of the parameter c (color palette).

The intersection points with section Σ 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 ) Γ u equals v s := ( x s , y s , w s ) Γ s. We compute the euclidean norm of the difference v u 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 4. Notice that this value coincides with the one obtained using the parameterization method and corresponds to a wave velocity of 6.14 mm/min. In Fig. 12(e), we display the heteroclinic connection between p l 1 and p r resulting for c ~.

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 v e l c 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 / d t = d V A / d t = 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 ε 0, 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 δ 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 ε 0 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 ε 0 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 theory12,13 to describe the behavior of the solutions close to the fold of the slow manifold and techniques from computer assisted proofs24 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 work10 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).

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/)

The authors have no conflicts to disclose.

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).

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

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

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 Δ x, for i = 1 , , N and Δ x is the distance between neurons.

The diffusion terms are approximated by central second-order finite differences, for instance, in the [ K + ] e case,
2 [ K + ] e x 2 [ K + ] e i + 1 2 [ K + ] e i + [ K + ] e i 1 Δ x 2 , i = 2 , , N 1 ,
(A1)
where i indicates the array and N denotes, therefore, its total number of elements (i.e., number of neurons and astrocytes) considered in the discrete network.
Dirichlet boundary conditions have been applied to the end points of the array, which aims to simulate healthy tissue. Thus, for the first and last indices, we must impose such conditions to obtain valid approximations there. Indeed,
First index approximation:  2 [ K + ] e x 2 [ K + ] e 2 2 [ K + ] e 1 + [ K + ] e , bdry Δ x 2 , Last index approximation:  2 [ K + ] e x 2 [ K + ] e , bdry 2 [ K + ] e N + [ K + ] e N 1 Δ x 2 .
(A2)
Same approximations apply to the diffusion term in the [ N a + ] e equation for system (1).

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

In this section, we present the computations of the derivatives of the functions H ( z ), for = l , m , r given in (22) and plotted in Fig. 6. Thus,
d H d z ( z ) = h x ( X ( z ) , Y ( z ) , z ) d X d z ( z ) + h y ( X ( z ) , Y ( z ) , z ) d Y d z ( z ) + h z ( X ( z ) , Y ( z ) , z ) ,
(B1)
where h ( x , y , z ) is the function given in (14). The first partial derivatives of h are
h x = 10 S N F Ω e ( g K 4 n 3 n ( x V K ) + g K n 4 ) , h y = 10 S A F Ω e P K F 2 R T ( z e ϕ [ K + ] i A e ϕ 1 + ϕ e ϕ z [ K + ] i A ( e ϕ 1 ) 2 ) , h z = 10 S N F Ω e ( g K n 4 R T F [ K + ] i + α 0 ( 1 + Ω a / Ω n ) z z [ K + ] i 4 ρ N z K K ( z + K K ) 3 ( [ N a + ] i K N a + [ N a + ] i ) 3 ) + 10 S A F Ω e ( P K F ϕ e ϕ e ϕ 1 4 ρ A z K K ( z + K K ) 3 ( [ N a + ] i A K N a , A + [ N a + ] i A ) 3 ) ,
(B2)
where ϕ = y F R T. For ϕ = 0, we use the L’Hôpital’s rule to make the derivatives continuous.
The expressions for the derivatives of the functions X ( z ) and Y ( z ) are obtained from the differentiation of the equations f ( X ( z ) , z ) = 0, for = l , m , r, and g ( Y ( z ) , z ) = 0 with respect to z. Indeed,
d X d z = f z f x and d Y d z = g z g y .

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

Developing functions f, g, and h as Taylor series around the coefficients of order 0, m 0 and n 0, yields
f ( m , z ) = f ( m 0 , z ) + f x ( m 0 , z ) ( ε m 1 + ε 2 m 2 ) + 1 2 2 f x 2 ( m 0 , z ) ε 2 m 1 2 + O ( ε 3 ) , g ( n , z ) = g ( n 0 , z ) + g y ( n 0 , z ) ( ε n 1 + ε 2 n 2 ) + 1 2 2 g y 2 ( n 0 , z ) ε 2 n 1 2 + O ( ε 3 ) , h ( m , n , z ) = h ( m 0 , n 0 , z ) + h x ( m 0 , n 0 , z ) ε m 1 + h y ( m 0 , n 0 , z ) ε n 1 + O ( ε 2 ) .
(C1)

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

  1. O(1):
    f(m0,z)=0m0(z,w)=Xr(z)(see the red curve in {Fig. 5}),g(n0,z)=0n0(z,w)=Y(z)(see the orange curve in {Fig. 5}).
  2. O(ε):
    Observe that, since
    f(Xr(z),z)=0
    and g(Y(z),z)=0, we have used
    m0z=dXr(z)dz=fz(m0,z)fx(m0,z)andn0z=dYdz=gz(n0,z)gy(n0,z).
    (C2)
  3. O(ε2) :
where
m1z=cw(fxz(m0,z)m0z+fzz(m0,z))fx(m0,z)2fz(m0,z)(fxx(m0,z)m0z+fzx(m0,z))(fx(m0,z))3,n1z=cw(gyz(n0,z)n0z+gzz(n0,z))gy(n0,z)2gz(n0,z)(gyy(n0,z)n0z+gzy(n0,z))(gy(n0,z))3.
(C3)
In this section, we detail the expression of the first partial derivatives of the functions f ~ ( x , z ) and g ~ ( y , z ) given in (14). The derivatives of higher order can be computed analogously. Recall that x = V N, y = V A, and z = [ K + ] e. The term 1 / C m has been omitted because, in practice, C m is equal to 1. Thus,
f ~ x = [ g N a 3 m 2 m ( 1 n ) ( x V N a ) + g N a m 3 ( 1 n x n + n V N a ) + g N a P m h p ( x V N a ) + g N a P m p h p + g K 4 n 3 n ( x V K ) + g K n 4 + g L ] , f ~ z = [ g K n 4 R T F [ K + ] i + α 0 ( 1 + Ω a / Ω n ) z z [ K + ] i + 2 ρ N z K K ( z + K K ) 3 ( [ N a + ] i K N a + [ N a + ] i ) 3 ] , g ~ y = [ P N a F 2 R T ( [ N a + ] e e ϕ [ N a + ] i A e ϕ 1 + ϕ e ϕ [ N a + ] e [ N a + ] i A ( e ϕ 1 ) 2 ) + P K F 2 R T ( z e ϕ [ K + ] i A e ϕ 1 + ϕ e ϕ z [ K + ] i A ( e ϕ 1 ) 2 ) ] , g ~ z = [ 2 ρ A z K K , A ( z + K K , A ) 3 ( [ N a + ] i A K N a , A + [ N a + ] i A ) 3 + P K F ϕ e ϕ e ϕ 1 ] ,
(D1)
where ϕ := y F R T.
We also define the derivatives when y = 0 using the L’Hôpital’s rule, and we obtain
g ~ y ( y , z ) | y = 0 = 1 2 F 2 R T ( P N a ( [ N a + ] e + [ N a + ] i A ) + P K ( z + [ K + ] i A ) )
and
g ~ z ( y , z ) | y = 0 = [ 2 ρ A z K K , A ( z + K K , A ) 3 ( [ N a + ] i A K N a , A + [ N a + ] i A ) 3 P K F ] .
To compute the eigenvalues and first-order approximation of the 1D stable manifold Γ ~ s, we need to compute the differential matrix for the restricted system (38), referred as F 2, given by
D F 2 ( z , w ) = ( 0 ε ε ( h x x z + h y y z + h z ) ε ( c ( h x x w + h y y w ) ) ) ,
(E1)
where x and y are defined as functions of z and w in (37). Thus,
x z = z ( m 0 + ε m 1 + ε 2 m 2 ) = f z ( m 0 , z ) f x ( m 0 , z ) + ε m 1 z + ε 2 m 2 z , x w = w ( m 0 + ε m 1 + ε 2 m 2 ) = ε c f z ( m 0 , z ) ( f x ( m 0 , z ) ) 2 + ε 2 m 2 w , y z = z ( n 0 + ε n 1 + ε 2 n 2 ) = g z ( n 0 , z ) g y ( n 0 , z ) + ε n 1 z + ε 2 n 2 z , y w = w ( n 0 + ε n 1 + ε 2 n 2 ) = ε c g z ( n 0 , z ) ( g y ( n 0 , z ) ) 2 + ε 2 n 2 w ,
where the derivatives m 0 / z and n 0 / z are taken from (C2) and we have used m 1 / w = n 1 / w = 0. Finally, m 1 / z and n 1 / z are given in (C3). For the other ones, we have
m 2 z = 1 2 [ ( f x x x ( m 0 , z ) z m 0 + f z x x ( m 0 , z ) ) f x ( m 0 , z ) f x x ( m 0 , z ) ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) ( f x ( m 0 , z ) 2 m 1 2 + f x x ( m 0 , z ) f x ( m 0 , z ) 2 m 1 m 1 z ] + [ c w ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) ( f x ( m 0 , z ) ) 2 m 1 z + c w f x ( m 0 , z ) 2 m 1 z 2 ] c [ ( h x ( m 0 , n 0 , z ) z m 0 h y ( m 0 , n 0 , z ) z n 0 h z ( m 0 , n 0 , z ) ) c f z ( m 0 , z ) ( f x ( m 0 , z ) ) 3 + ( c w h ( m 0 , n 0 , z ) ) c ( f x z ( m 0 , z ) z m 0 + f z z ( m 0 , z ) ) f x ( m 0 , z ) 3 f z ( m 0 , z ) ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) ( f x ( m 0 , z ) ) 4 ] , m 2 w = 1 2 f x x ( m 0 , z ) f x ( m 0 , z ) 2 m 1 ( c f z ( m 0 , z ) ( f x ( m 0 , z ) ) 2 ) + 2 c f x ( m 0 , z ) m 1 z c 3 f z ( m 0 , z ) ( f x ( m 0 , z ) ) 3 ,
with
2 m 1 z 2 = ( c w ) 1 ( f x ( m 0 , z ) ) 4 { [ ( ( f x x z ( m 0 , z ) z m 0 + f z x z ( m 0 , z ) ) z m 0 + f x z ( m 0 , z ) z z 2 m 0 + ( f x z z ( m 0 , z ) z m 0 + f z z z ( m 0 , z ) ) ) f x ( m 0 , z ) ( f x z ( m 0 , z ) z m 0 + f z z ( m 0 , z ) ) ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) 2 f z ( m 0 , z ) ( ( f x x x ( m 0 , z ) z m 0 + f z x x ( m 0 , z ) ) z m 0 + f x x ( m 0 , z ) z z 2 m 0 + ( f x z x ( m 0 , z ) z m 0 + f z z x ( m 0 , z ) ) ) ] f x ( m 0 , z ) 3 ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) ( ( f x z ( m 0 , z ) z m 0 + f z z ( m 0 , z ) ) f x ( m 0 , z ) 2 f z ( m 0 , z ) ( f x x ( m 0 , z ) z m 0 + f z x ( m 0 , z ) ) ) } .
One obtains the expressions for the partial derivatives of n 2 with respect to z and w by carefully exchanging f’s and x’s for g’s and y’s, respectively.
1.
J. P.
Dreier
, “
The role of spreading depression, spreading depolarization and spreading ischemia in neurological disease
,”
Nat. Med.
17
,
439
447
(
2011
).
2.
D. R.
Kramer
,
T.
Fujii
,
I.
Ohiorhenuan
, and
C. Y.
Liu
, “
Cortical spreading depolarization: Pathophysiology, implications, and future directions
,”
J. Clin. Neurosci.
24
(
3
),
22
27
(
2016
).
3.
D.
Pietrobon
and
M. A.
Moskowitz
, “
Chaos and commotion in the wake of cortical spreading depression and spreading depolarizations
,”
Nat. Rev. Neurosci.
15
(
6
),
379
393
(
2014
).
4.
M.
Lauritzen
,
J. P.
Dreier
,
M.
Fabricius
,
J. A.
Hartings
,
R.
Graf
, and
A. J.
Strong
, “
Clinical relevance of cortical spreading depression in neurological disorders: Migraine, malignant stroke, subarachnoid and intracranial hemorrhage, and traumatic brain injury
,”
J. Cereb. Blood Flow Metab.
31
(
1
),
17
35
(
2011
).
5.
G. G.
Somjen
, “
Mechanisms of spreading depression and hypoxic spreading depression-like depolarization
,”
Physiol. Rev.
81
(
3
),
1065
1096
(
2001
).
6.
A.
Bellot-Saez
,
O.
Kkesi
,
J. W.
Morley
, and
Y.
Buskila
, “
Astrocytic modulation of neuronal excitability through K + spatial buffering
,”
Neurosci. Biobehav. Rev.
77
,
87
97
(
2017
).
7.
B.
Ma
,
R.
Bucklew
,
Y.
Du
,
C. M.
Kiyoshi
,
C. C.
Alford
,
W.
Wang
,
D. D.
McTigue
,
J. J.
Enyeart
,
D.
Terman
, and
M.
Zhou
, “
Gap junction coupling confers isopotentiality on astrocyte syncytium
,”
Glia
64
(
2
),
214
226
(
2016
).
8.
M.
Nedergaard
and
U.
Dirnagl
, “
Role of glial cells in cerebral ischemia
,”
Glia
50
,
281
286
(
2005
).
9.
Y.
Zhao
and
D. A.
Rempe
, “
Targeting astrocytes for stroke therapy
,”
Neurotherapeutics
7
,
439
451
(
2010
).
10.
G.
Huguet
,
A.
Joglekar
,
L. M.
Messi
,
R.
Buckalew
,
S.
Wong
, and
D.
Terman
, “
Neuroprotective role of gap junctions in a neuron astrocyte network model
,”
Biophys. J.
111
(
2
),
452
462
(
2016
).
11.
R.
Lee
, “Analysis of spreading depolarization as a traveling wave in a neuron-astrocyte network,” M.Sc. thesis (The Ohio State University, 2017).
12.
C. K. R. T.
Jones
, “Geometric singular perturbation theory,” in Dynamical Systems. Lecture Notes in Mathematics, edited by R. Johnson (Springer, Berlin, 1995), Vol. 1609.
13.
C.
Kuehn
, “Geometric singular perturbation theory,” in Multiple Time Scale Dynamics. Applied Mathematical Sciences (Springer, Cham, 2015), Vol. 191.
14.
X.
Cabré
,
E.
Fontich
, and
R.
De La Llave
, “
The parameterization method for invariant manifolds III: Overview and applications
,”
J. Differ. Equ.
218
(
2
),
444
515
(
2005
).
15.
A.
Haro
,
M.
Canadell
,
J. L.
Figueras
,
A.
Luque
, and
J. M.
Mondelo
,
The Parameterization Method for Invariant Manifolds
(
Springer
,
Berlin
,
2016
).
16.
N.
Fenichel
, “
Persistence and smoothness of invariant manifolds for flows
,”
Indiana Univ. Math. J.
21
,
193
226
(
1971/1972
).
17.
N.
Fenichel
, “
Asymptotic stability with rate conditions
,”
Indiana Univ. Math. J.
23
,
1109
1137
(
1973/1974
).
18.
M.
Hirsch
,
C.
Pugh
, and
M.
Shub
, “Invariant manifolds,” in Lecture Notes in Mathematics (Springer, Berlin, 1977), Vol. 538.
19.
B. G.
Ermentrout
and
D. H.
Terman
,
Mathematical Foundations of Neuroscience
(
Springer
,
New York
,
2010
).
20.
K.
Keener
and
J.
Sneyd
, “Mathematical physiology,” in Interdisciplinary Applied Mathematics (Springer, 2008).
21.
A.
Griewank
and
A.
Walther
,
Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation
(
SIAM
,
2008
), Vol. 105.
22.
G.
Huguet
and
R.
de la Llave
, “
Computation of limit cycles and their isochrons: Fast algorithms and their convergence
,”
SIAM J. Appl. Dyn. Syst.
12
(
4
),
1763
1802
(
2013
).
23.
A.
Pérez-Cervera
,
T.
M-Seara
, and
G.
Huguet
, “
Global phase-amplitude description of oscillatory dynamics via the parameterization method
,”
Chaos
30
(
8
),
083117
(
2020
).
24.
M. J.
Capiński
,
J.
Gonzalez
,
J. P.
Marco
, and
J. D.
Mireles James
, “
Computer assisted proof of drift orbits along normally hyperbolic manifolds
,”
Commun. Nonlinear Sci. Numer. Simul.
106
,
105970
(
2022
).
25.
H. P.
McKean
, “
Nagumo’s equation
,”
Adv. Math.
4
,
209
223
(
1970
).
26.
J.
Rinzel
and
J. B.
Keller
, “
Traveling wave solutions of a nerve conduction equation
,”
Biophys. J.
13
,
1313
1337
(
1973
).
27.
L.
Lemaire
,
M.
Desroches
,
M.
Krupa
, and
M.
Mantegazza
, “
Idealized multiple-timescale model of cortical spreading depolarization initiation and pre-epileptic hyperexcitability caused by NaV1.1/SCN1A mutations
,”
J. Math. Biol.
86
(
6
),
92
(
2023
).
Published open access through an agreement with Universitat Politècnica de Catalunya