Traveling waves in a model for cortical spreading depolarization with slow-fast dynamics

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 to 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. , 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 an 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.


Introduction
Cortical spreading depression and depolarization (CSD) are slowly propagating waves of rapid, nearcomplete depolarization of neurons and glial cells across the cortex [4,24,15].These waves move slowly through the cortex (2 to 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 [17,15].
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 [26,4,15].Astrocytes may play a crucial role in regulating potassium concentration by absorbing its accumulation from the extracellular space [1], thus potentially delaying or even preventing the onset of CSD waves [21,27,22].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 occur 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 [12].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 [12] 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 [12] based on the methodology developed in [18].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.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 [16,13] and introduce a small parameter that represents the ratio between the different time scales.Using techniques from classical singular perturbation theory [16,13], we break down the system into its fast and slow subsystems and construct an 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 [2,9] and another employing Fenichel's theory, which ensures the persistence of normally hyperbolic invariant manifolds [6,7,10].
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 Section 2 we present the original model introduced in [12] to describe CSD and the reduced version capturing the initiation of the wave.In Section 3 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 Section 4 we provide the computation of the heteroclinic orbit and the velocity for the reduced system using the parameterization method [9,2].In Section 5 we provide an alternative method to compute the heteroclinic orbit numerically based on Fenichel's theory [6,7].We end with a discussion in Section 6.Finally, the Appendix contains some technical details involving algebraic mathematical computations.
2 The mathematical model

Model for cortical spreading depression involving astrocytes
In this section, we present the neuron-astrocyte network model introduced in [12], 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 (N a + ) 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 [4].A thorough study on preventing or delaying this phenomenon was carried out in [12], 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, e , respectively, and the gating variables n and h p .The system of differential equations is given by: The ionic currents I N a , I N aP , I K and I L are modelled as in the Hodgkin-Huxley model, that is, (2) The reversal potential E X , for X ∈ {N a, K}, is modelled by the Nernst equation where R, T , and F are the gas constant, temperature, and Faraday's constant, respectively.The values are given in Table 2.
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 for Finally, the time constants modeling the rate at which the gating variables reach its steady state values are The ionic currents I A N a and I A K in the modelization of the astrocyte equation are described by the Goldman-Hodgkin-Katz (GHK) formulation: with φ := F RT V A , P N a = 1.5 • 10 −8 and P K = 1 • 10 −6 .Notice that the currents I A N a and I A K 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 A N a and I A K at V A = 0, 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 modelled using the GHK formalism in [12].Since in this paper we will not consider them, we do not include details of the modelization here and we refer the reader to [12] for further information.
The sodium-potassium pump currents in neurons and astrocytes, I P m and I A P m , respectively, are given by where ρ N , ρ A represent the maximal current allowed to go through the neuron and astrocyte pumps, respectively (regarded as parameters in [12] 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 [12] that a wave can be initiated and propagated through the tissue.We know from previous works [12,18] 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 2 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 endpoints 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 −30mV.The simulation shows a depolarization wave which propagates throughout the array of neuron-astrocyte pairs (see Figure 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 ± ct), whose spatial profile advances in time at a constant speed c (therefore being invariant by temporal translations).

Model reduction
In this section, we present a reduction of model (1) following the strategy employed in [18].First, we provide details of the simplifications and its main differences with respect to [18].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 + ] A i 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 + ] A i = 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 + ] A i .Thus, the variable [K + ] A i in equation ( 1) is also taken as constant.Finally, the conservation of Figure 1: Solutions of the model (1) with simulated K + injection.The simulation consists of a 1D array network of 50 neuron-astrocyte pairs sharing the same extracellular space, with the middle four cells externally stimulated with an insult of K + at a rate of 5 mM/s until a depolarizing wave is triggered.The initial conditions are identical for all the pairs.The diffusive terms in system (1) are approximated by central second-order finite differences, imposing Dirichlet boundary conditions on the endpoints of the array (using typical levels of K + and N a + in healthy tissues).(Top) Temporal profile of the variable V N for each cell during a CSD wave (each curve is coloured according to the neuron index).(Bottom) 2D views of the traveling wave solution (arising) on V N , V A , [K + ] e and [N a + ] e , spreading through space from middle cells towards adjacent ones over time.Parameters are given in Table 2.
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 + ] A i 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 [12,21] 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 Figure 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: System (8) is the model we will analyze in this paper, which is the same considered in [18] except that we keep the terms I A K and I A P m in the equation for [K + ] e , describing the influence of astrocytes in the extracellular levels of K + .
In Figure 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 Figure 1) but traveling fronts.As in simulations of Figure 1, we modelled 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 Figure 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 1).
A further simplification of model ( 8) has been done in [18], 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 Figure 2).Mathematically, imposing that the system is reduced to a single PDE for [K + ] e of reaction-diffusion type, that preserves qualitatively the dynamics shown in Figure (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.
3 Traveling wave solutions in the singular limit 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  [12] and [18].Note that E X in (3) must be multiplied by a factor 1000 to convert from V to mV.
where the constant c (to be determined) establishes the wave propagation velocity, vel = c √ D K .Substituting (10) in ( 8) results in a system of ordinary differential equations, given by We will see that system (11) has two saddle points p l1 and p r (see Table 3).To find a traveling front solution of system (8), we will look for an heteroclinic orbit of system (11) from p l1 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 := = dz dξ , to finally obtain which, for simplicity, can be written in a more general form with where we have explicitly indicated the dependence on x, y and z for each ionic current.for fixed values of y (as indicated in each panel).The zero-level curve for each function is displayed as the intersection of the grey plane with each surface.Notice also the different orders of magnitude on the value ranges that each function takes, suggesting distinct scales present in system (16).Parameters of the functions are taken from Table 2.
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 Figure 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 dx dξ and dy 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, where we have denoted f = ε f and g = ε g.Equivalently, re-scaling the time variable εη = ξ, system (15) becomes Remark 3.1.By looking at the range of values taken by g in Figure 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 ẋ (to indicate the dynamics of x are 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 equations ( 15) and ( 16) correspond to a singular perturbation problem in parameter ε.First, we are going to construct an heteroclinic orbit in the singular limit, that is when ε → 0 by splitting the dynamics into the slow and fast dynamics.

Dynamics of the slow subsystem
Taking ε = 0 in system (15), we obtain the so-called slow subsystem -a system of differentialalgebraic equations -given by 0 The first two algebraic equations (17a)-(17b) describe the commonly named critical manifold M 0 , a 2-dimensional manifold defined as where Figure 5A 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, 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 Figure 5A).Thus, the manifold M 0 has three branches (see Figure 5B), namely as well as two fold lines In Figure 5B 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 grey curve) and a solution of the simulation of system (8) portrayed in Figure 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 l 0 and the upper branch M r 0 , 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 Section 3.3.
Figure 5: Dissection of the critical manifold defined by (17a)-(17b).(A) Zero curves of f (x, z) (red) and g(y, z) (orange) on the (x, z) and (y, z) planes, respectively.(B) The manifold M 0 (dashed grey curve) defined in (19), which is the projection of the 2D critical manifold M 0 onto the (z, x, y) space, and a representative solution of the simulation of system (8) displayed in Figure 2 (cyan curve starting at the black-filled dot).
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 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 for * = l, m, r.Notice that each function H * is defined for a different range of z.In Figure 6 we plot the functions H * (z) (distinguished by the colour), each one defined in its respective domain.
It turns out that H l has two zeros, z l 1 and z l 2 , 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 3.Two of the equilibrium points, p l1 = (X l (z l 1 ), Y (z l 1 ), z l 1 , 0) and p l2 = (X l (z l 2 ), Y (z l 2 ), z l 2 , 0), lie on the lower branch M l 0 of the manifold and the third one, p r = (X r (z r ), Y (z r ), z r , 0), on the upper branch M r 0 .
If we denote by F * ss (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  To compute them we have applied a Newton method to solve the system (f, g, h)(x, y, z) = 0.The expressions for the first partial derivatives are given in (46) and (42).
Observe that det DF * ss = dH * dz and Tr DF * ss = c > 0. The expressions for the derivatives of H * are given in Appendix B. Function H l has negative slope at point z l 1 and so does H r at z r , hence they are saddle points.The derivative of H l at z l 2 is positive, and thereby assuring its unstable character.To illustrate this, in Figure 7 we plot both eigenvalues λ ± for each equilibrium point for a range of positive values c.

Dynamics of the fast subsystem
Setting now ε = 0 in system (16), we obtain the so-called 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) ∈ R 4 associated to 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 Figure 8 we plot these eigenvalues as functions of z.
for each branch X * (z), for * = l, m, r : λ l 1 (blue), λ m 1 (purple) and λ r 1 (red), while on the right axis, we show in orange the second eigenvalue, λ 2 (z) = 1 c g y (Y (z), z).Vertical dashed lines indicate the z-coordinate of the folds of the critical manifold (see Figure 5).Inset shows a close-up of the region comprised by the folds, where the three solutions X * (z), for * = l, m, r coexist.
Observe that λ 2 < 0 for all values of z.For points on M l 0 , the eigenvalue λ l 1 < 0 and therefore, both eigenvalues are negative.The same happens for points on M r 0 .However, for points on M m 0 , λ m 1 > 0. From this analysis we conclude that the lower and upper branches M l 0 and M r 0 are normally hyperbolic attracting manifolds, while the middle branch M m 0 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.

Singular heteroclinic connection concatenating slow and fast subsystems
In this section we look for an heteroclinic connection between points p l1 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 l 0 until it reaches the fold F R , where it rapidly switches to track the upper branch M r 0 , we build a piecewise system defined on the (z, w) plane of the form where H(z) is a piecewise smooth function defined as where z R = 18.276 is the z-coordinate of the right fold F R .We look for an heteroclinic orbit between the two saddle points, p l1 = (z l 1 , 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 l1 and p r , respectively, which are 1D curves, and extend them until reaching the section where the definition of function H(z) changes from the lower branch M l 0 to the upper branch M r 0 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 Figure 9 we plot, for different values of c, the unstable and stable invariant manifolds of p l1 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.07426.
We recognize (24) as the bistable reaction-diffusion equation studied in [18].Adapting classical methods described in [14,5], one can show that the function H(z) satisfies the hypothesis that guarantee the existence of heteroclinic orbits for system (24).4 Computing the heteroclinc orbit for system (15) with ε = 0 by means of the parameterization method 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 r 0 and M l 0 .Fenichel's theory [6,7] 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 l1 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 l1 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 l1 , 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 The singular heteroclinic orbit (black curve) happens for c 0 ≈ 0.07426.Notice the different ranges of the left and right intervals.This has been done for visualization purposes, so that the first part of the heteroclinic orbit can be fully appreciated.
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 an 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 3dimensional, 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 an 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 l1 Γ 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 method [9] 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 Simultaneously, we will also determine the internal dynamics of the manifold, governed by the ODE ṡ = f (s) with f : R → R. Imposing W (s) to be a solution of the system Ẋ = F (X) (system (15) in our case), we obtain the invariance equation Assume W (s) and f (s) can be both expressed as power series in s, 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 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 (29) The equation of order 1 leads to the eigenvalue problem thus begin 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 which corresponds to With this simplification, the equation of order 2 reads as where denotes the coefficient of order 2 of the Taylor expansion in s.The resulting system is linear 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 where, again, the coefficient F k has been split into the known, F (W <k ) | k , and unknown parts, DF (p r )W k .
Equation (34) can be solved provided det(DF (p r ) − kλ slow Id) = 0 , or, equivalently, that kλ slow is not an eigenvalue of the differential matrix DF (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 equation ( 34) can be numerically computed efficiently with accuracy as high as desired using the Automatic Differentiation algorithms [9,8].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 [9]).27) as a function of the parameter s for a parametrization W (s) computed up to order k = 55 with the internal dynamics f (s) given by (31).Any value s whose error is smaller than the allowable threshold (i.e.below the dashed-black line) produces a potentially accurate initial condition to be backward-integrated.
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)) − DW (s)f (s)|| ∞ , is less than 10 −10 (see Figure 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 l1 , Γ u we integrate the system forward starting from a point obtained using a first order approximation of the manifold (i.e.p l1 + δ v l1 with v l1 the eigenvector of the linearized system around p l1 associated to the positive eigenvalue and δ a small value).
In Figure 11A-D we illustrate the strategy used for computing an heteroclinic orbit.For a range of c-values including the value c 0 = 0.07426 (which is the value for which we found an 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 Figure 11A and 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 Figure 11B and D).For an 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 blue-dotted line in Figure 11F).In Figure 11F 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 vel = c .96 • 0.073135 mm/s = 0.102389 mm/s = 6.1433 mm/min.

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 Γ 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 r 0 .By doing so, we assure 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 l1 , Γ 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 r 0 , it can be approximated by considering both fast variables x and y as formal expansions in ε, 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 equation ( 20)).
From imposing the manifold (35) to be invariant by system (15) we reach the following system of equations: 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: and and Next, we will restrict the dynamics of the full system (15) onto the second order approximation of the manifold M r ε by taking The resulting system is ż = εw , 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 Figure 12A and 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 grey), 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 Figure 12B and 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 Figure 12F).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.073135 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 Figure 12E we display the heteroclinic connection between p l1 and p r resulting for c.

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 [12], 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 [13,16], and numerically efficient methods based on the parameterization method [9], 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 to the presence of different time-scales in the system.
We emphasize that the computation of the propagation velocity vel ≈ c √ D K by means of computing an 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 [18].However, we want to highlight some relevant differences with respect to this work.In [18] the system is further reduced to a single PDE for [K + ] e of reaction-diffusion type, imposing that the dynamics for V N and V A in (8) are set to be instantaneous (i.e.dV n /dt = dV 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 ε = 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 Figures 11AC and 12AC).Indeed, in Figure 13 we present the numerical simulations of (8) taking V N and V A as instantaneous variables (i.e.dV N dt = dV A dt = 0).As expected from our results, these simulations reveal a good agreement with those for system (8) shown in Figure 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 Figure 4 (see Remark 3.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 [2,9], widely used in the field of dynamical systems [11,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 theory [13,16] to describe the behaviour of the solutions close to the fold of the slow manifold and techniques from computer assisted proofs [3] 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 [12] 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 [12], the increase in [K + ] e is accompanied by a steep decline in [N a + ] e (see Figure 1).Moreover, freezing the inactivation variable h p of the persistent sodium current I N aP 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 [25,20,14].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 [19] incorporating the spatial domain).
Figure 13: Simulations as in Figure 2 for the reduced model ( 8) with V N and V A now regarded as variables with instantaneous dynamics, that is, setting dV N dt = dV A dt = 0 on system (8).Variables V N and V A are described in terms of the extracellular concentration of K + by means of the algebraic constraints given in (9).Comparing with the simulation in Figure 2, the dynamics of the model ( 8) is close to that of the limit case considered here.

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

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) (Figures 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 where i indexes the array and N denotes therefore its total number of elements (i.e.number of neurons and astrocytes) considered in the discrete network.Dirichlet boundaries conditions have been applied to the endpoints of the array, which aims to simulate healthy tissue.Thus, for the first and last indexes, we must impose such conditions to obtain valid approximations there.Indeed, First index approximation: Last index approximation: 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.

B Computation of the derivatives of the functions H *
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 Figure 6.Thus, where h(x, y, z) is the function given in (14).The first partial derivatives of h are where φ = yF RT .For φ = 0 we use 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, dX * dz = − ∂f /∂z ∂f /∂x , and dY dz = − ∂g /∂z ∂g /∂y .

C Taylor expansions of the manifold M r ε in ε
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 After substituting the Taylor expansions (43) into system (36), we are able to solve the resulting system of equations up to any power of ε: • O(1) : f (m 0 , z) = 0 =⇒ m 0 (z, w) = X r (z) , (see red curve in Figure 5) Observe that, since f (X r (z), z) = 0 and g(Y (z), z) = 0, we have used where ∂m 1 ∂z = −cw (f xz (m 0 , z) ∂m0 ∂z + f zz (m 0 , z)) f x (m 0 , z) − 2f z (m 0 , z)(f xx (m 0 , z) ∂m0 ∂z + f zx (m 0 , z)) (f x (m 0 , z)) 3, ∂n 1 ∂z = −cw (g yz (n 0 , z) ∂n0 ∂z + g zz (n 0 , z)) g y (n 0 , z) − 2g z (n 0 , z)(g yy (n 0 , z) ∂n0 ∂z + g zy (n 0 , z)) (g y (n 0 , z)) 3  , D Computation of the partial derivatives of the functions f (x, z) and g(y, z) 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, where φ := yF RT .We also define the derivatives when y = 0 using L'Hôpital's rule, and we obtain 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 where x and y are defined as functions of z and w in (37).Thus, where the derivatives ∂m 0 /∂z and ∂n 0 /∂z are taken from (44) and we have used ∂m 1 /∂w = ∂n 1 /∂w = 0. Finally, ∂m 1 /∂z and ∂n 1 /∂z are given in (45).For the other ones, we have 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.

Figure 2 :
Figure 2: Solutions of model (8) with simulated K + injection.Fifty neuron-astrocyte pairs are simulated, with the middle four cells stimulated with insult of K + at a rate of 5 mM/s until a front of spreading depolarization is evoked.(Top) Temporal profile of the variables V N , V A and [K + ] e for each cell during a front of spreading depolarization (each curve is coloured according to the neuron-astrocyte index).(Bottom) 2D views of a depolarizing front on V N , V A and [K + ] e spreading through space to each neighbouring cell over time.As in Figure1, Dirichlet boundary conditions have been applied on the extracellular potassium K + , to simulate typical levels of [K + ] e in healthy tissues.Parameters are given in Table2.

20 Figure 3 :
Figure3: Estimation of the wave speed of the numerical simulations of system(8) shown in Figure2from the voltage traces of neurons 10 and 20.Using that the distance between cells is ∆x = 0.044 mm, the estimated front velocity is vel ≈ 3.82 mm/min.

Figure 4 :
Figure 4: Plots of the graphs of functions given in (14): (A) f (x, z), (B) g(y, z) and (C-F) h(x, y, z)for fixed values of y (as indicated in each panel).The zero-level curve for each function is displayed as the intersection of the grey plane with each surface.Notice also the different orders of magnitude on the value ranges that each function takes, suggesting distinct scales present in system(16).Parameters of the functions are taken from Table2.

Figure 6 :
Figure 6: Graph of functions H * (z) := h X * (z), Y (z), z for * = l (blue), m (purple) and r (red).The filled (resp.empty) dots indicate the zeros of the functions H * corresponding to stable (resp.unstable) points.Inset shows a zoom of the region close to the zeros z l1 and z l 2 (lying very close to each other).To compute them we have applied a Newton method to solve the system (f, g, h)(x, y, z) = 0.The expressions for the first partial derivatives are given in (46) and (42).

Figure 7 :
Figure7: Eigenvalues of the slow subsystem(17) associated to each of its equilibrium points as c varies from 0 to 0.1.Each color corresponds to an equilibrium point and the eigenvalues are differentiated by the line-style: solid lines for λ + and dashed ones for λ − .Eigenvalues for p l2 are complex conjugates for small values of c and we display only the real part of them.For the range of c shown, points p l1 and p r are of saddle-type while p l2 is unstable.

Figure 8 :
Figure 8: Eigenvalues for the equilibrium points of the fast subsystem (23a)-(23b) in terms of (parameter) z for a value c 1. On the left y-axis we depict the eigenvalue λ * 1(z) = 1 c f x (X * (z), z)for each branch X * (z), for * = l, m, r : λ l 1 (blue), λ m 1 (purple) and λ r 1 (red), while on the right axis, we show in orange the second eigenvalue, λ 2 (z) = 1 c g y (Y (z), z).Vertical dashed lines indicate the z-coordinate of the folds of the critical manifold (see Figure5).Inset shows a close-up of the region comprised by the folds, where the three solutions X * (z), for * = l, m, r coexist.

Figure 9 :
Figure9: Computation of the heteroclinic orbit in the (singular) slow subsystem(24).Blue curves depict the unstable invariant manifolds of p l1 for two different values of c up to the Poincaré section Σ, located at fold z R = 18.276 (dashed-grey line).Orange curves represent the stable invariant manifolds of p r for the same two values of c.The curves are distinguished by the degree of transparency; the less transparent line corresponds to c = 0.04 while the opaque one to c = 0.09.The singular heteroclinic orbit (black curve) happens for c 0 ≈ 0.07426.Notice the different ranges of the left and right intervals.This has been done for visualization purposes, so that the first part of the heteroclinic orbit can be fully appreciated.

Figure 10 :
Figure 10: Error of the invariance equation (27) as a function of the parameter s for a parametrization W (s) computed up to order k = 55 with the internal dynamics f (s) given by (31).Any value s whose error is smaller than the allowable threshold (i.e.below the dashed-black line) produces a potentially accurate initial condition to be backward-integrated.

Figure 11 :
Figure 11: Computation of an heteroclinic orbit between the two saddle points p l 1 and pr in the full system (15) using the parameterization method.(A,C) 3D projections of the phase space onto the: (A) (x, z, w) and (C) (y, z, w) coordinated systems.In light grey we depict the critical manifold M0 and in dark grey we represent the Poincaré section Σ1 = {z = 22}.The equilibrium points are denoted by circles, of which the black-filled ones correspond to the saddle points p l 1 and pr.The coloured orbits correspond to the unstable Γ u and (weakest) stable Γ s 1D manifolds of the saddle points for different values of c (colour palette).The 3D representation of the singular heteroclinic orbit, shown in Figure 9, has been represented by a blackdashed curve.(B,D) 2D projections of intersection points of the 1D manifolds Γ u (crosses) and Γ s (circles) with the Poincaré section Σ1 for the range of values c considered.Projection onto the: (B) (w, x) and (D) (w, y) planes.(E) 3D projection (onto the (z, w, x)-space) of the numerically computed heteroclinic connection for c. (F) Metrics computed upon the Poincaré section between the intersection points of the unstable and stable manifolds as a function of the parameter c: Euclidean distance (blue-dotted curve) and component-wise distances (see panel legend).

Figure 12 :
Figure 12: Computation of the heteroclinic orbit between the two saddle points p l 1 and pr when full system (15) is restricted to the Fenichel's manifold (35) on the upper branch.(A,C) 3D projections of the phase space onto the: (A) (z, w, x) and (C) (z, w, y) coordinated systems.The critical manifold M0 and the Poincaré section Σ1 = {z = 22} are depicted in light and dark gray colours, respectively.Black-filled (resp.empty) dots indicate the saddle (resp.unstable) equilibrium points of the system.The coloured curves correspond to the unstable Γ u and (weakest) stable Γ s 1D manifolds of the saddle points (up to the Poincaré section) for different values c ranging from 0.06 (blue) to 0.1 (green).The black-dashed curve illustrates the 3D representation of the singular heteroclinic orbit shown in Figure 9 (B,D) 2D projections of the intersection points of the 1D manifolds Γ u (crosses) and Γ s (circles) with the Poincaré section Σ1 for the range of values c considered.Projection onto the: (B) (w, x) and (D) (w, y) planes.(E) 3D projection (onto (z, w, x)-space) of the numerically computed heteroclinic connection for c. (F) Metrics computed upon the Poincaré section between the intersection points of the unstable, vu, and the stable, vs, manifolds as a function of the parameter c: Euclidean distance (blue-dotted curve) and component-wise distances (see panel legend).

=
[N a + ] e + [N a + ] A i + P K z + [K + ] − 2ρ A zK K,A (z + K K,A ) 3 [N a + ] A i K N a,A + [N a + ] A i E Linear local approximation for the 1D stable manifold Γ s F 2 , given by DF 2 (z, w)

Table 1 :
(8)imated wave speed of numerical simulations of system(8)as the number of neuronastrocyte pairs in the medium increases.

Table 2 :
Parameters for system (1) used in this paper which combine parameter values from N 5µA/cm 2 maximum neuron Na + -K + ATPase pump current ρ A 5 µA/cm 2 maximum astrocyte Na + -K + ATPase pump current