Nucleation of reaction-diffusion waves on curved surfaces

We study reaction-diffusion waves on curved two-dimensional surfaces, and determine the influence of curvature upon the nucleation and propagation of spatially localized waves in an excitable medium modelled by the generic FitzHugh-Nagumo model. We show that the stability of propagating wave segments depends crucially on the curvature of the surface. As they propagate, they may shrink to the uniform steady state, or expand, depending on whether they are smaller or larger, respectively, than a critical nucleus. This critical nucleus for wave propagation is modified by the curvature acting like an effective space-dependent local spatial coupling, similar to diffusion, thus extending the regime of propagating excitation waves beyond the excitation threshold of flat surfaces. In particular, a negative gradient of Gaussian curvature $\Gamma$, as on the outside of a torus surface (positive $\Gamma$), when the wave segment symmetrically extends into the inside (negative $\Gamma$), allows for stable propagation of localized wave segments remaining unchanged in size and shape, or oscillating periodically in size.


Introduction
Wave propagation in excitable extended media described by nonlinear reactiondiffusion equations has widespread applications in chemistry, biology, and medicine. A particularly important example are neuronal systems where spreading depression (SD) waves are associated with a pathological dysfunction of brain activity that occurs, for instance, during migraine or stroke. Understanding the effect of internal control mechanisms in neuronal wave dynamics is of great relevance not only for comprising the functionality of the human brain [1], but also for developing novel future therapies of pathological states which are connected with cortical spreading depression waves in migraine aura or stroke [2,3]. There is clinical and experimental evidence [4,5] that spatially localized wave segments play a dominant role in these phenomena. In general, spatially localized wave segments, as they propagate, might shrink, expand, or remain unchanged in size and shape, in which case they are called particle-like waves or dissipative solitons. Spatially localized wave segments also represent critical structures which can be stabilized by global feedback [6,7,8]. They play an important role for the nucleation of propagating waves and wave segments in two-dimensional spatial domains. Generally, waves can be controlled by feedback or closed loop control. This is a robust and versatile concept which uses the internal dynamics of the system to generate a control signal which directs the system towards a desired dynamics. A plethora of examples are provided by global or nonlocal and in some cases time-delayed feedback control of wave propagation in reaction-diffusion systems [9,10,11,12,13,14,15]. On the other hand, the curvature of the medium itself also provides a means of internal control of the stability, as we will show in this paper.
Most previous studies have focussed on wave propagation in planar spatial domains, yet there is also a considerable body of work on reaction-diffusion waves in curved surfaces mostly on spirals and ring waves [16,17,18,19,20,21,22,23,24] but not to the best of our knowlege on nucleation. The cortex, however, represents a strongly curved surface. It is the purpose of this paper to study nucleation and propagation of wave segments on curved two-dimensional surfaces. We will demonstrate that positive or negative Gaussian curvature of the spatial domain has a dramatically different effect upon the wave dynamics.
The paper is organized as follows. In Sect. 2 we present the model. In Sect. 3 we discuss wave solutions on a torus, which represents a curved surface on which locally both positive and negative Gaussian curvature occurs. We consider ring waves, wave segments, and dissipative solitons (critical nuclei) stabilized by feedback control. Specifically, we study ring wave break-up, curvature-induced changes of stability, and curvature-induced stabilisation of wave-segments. In Sect. 4 we draw conclusions. u v canard trajectory sub-threshold super-threshold Figure 1. The nullclinesu = 0 (solid red) andv = 0 (solid green) in the phase space of the homogeneous FHN system with β = 1.4. Their intersection at (u s , v s ) is a stable fixed point. Three trajectories are drawn for ε = 0.04: one canard-like trajectory (dotted), passing through the maximum of the nullclineu = 0, and two trajectories starting at v = v s nearby but on opposite sides of the canard trajectory. They diverge sharply, producing threshold behavior: the dashed and the dash-dotted trajectories represent super-threshold and sub-threshold stimulation, respectively.

Model
In this work we use the FitzHugh-Nagumo (FHN) system as a generic model for excitable systems [25], and study it in two spatial dimensions. It is a generalization of the van der Pol oscillator [26] and also known as Bonhoeffer-van der Pol oscillator [27,28]. The model was first suggested by FitzHugh in 1961 [29], and independently by Nagumo et al. in 1962 [30]. As a two-variable simplification of the four-variable Hodgkin-Huxley model, which describes the propagation of action potentials along the giant squid axon [31], it describes the response of an excitable nerve membrane to external current stimulation: The original interpretation of the FHN Eqs. (1)-( 2) is based on a single neuron: The variable u models fast changes of the electrical potential across the membrane of a nerve cell axon (occurring as spikes in the time series), and v is the recovery variable related to the gating mechanism of the membrane channels [29]. The small parameter ε 1 represents the time scale ratio of the two variables. The generic dynamical mechanism described by Eqs. (1)-( 2) was also suggested to describe a fundamentally different ionic Schematic phase diagram of different regimes of the FithHugh-Nagumo model in the plane of wavesize S and threshold β for fixed ε: weakly excitable (perturbations grow to spiral waves); subexcitable (perturbations shrink in length); nonexcitable (perturbations shrink in width, no propagation in one dimension). The respective boundaries are marked by ∂R ∞ and ∂P 1D . ∂R denotes the boundary of the critical nucleus of size S below which perturbations shrink. The red solid line marks the control loop, which stabilizes the critical nucleus (S * , β * ) indicated by a red dot. excitability of neuronal tissue that originates from bistable ion homeostasis and is the bais of spreading depression [32,33].
In general, the fast variable u is called the activator variable, whereas the slow variable v is referred to as inhibitor variable. The diffusion constant of the activator is D, which is chosen as 0.12 throughout this paper (and can simply be interpreted as a scaling of space), and inhibitor diffusion is assumed to be slow, and hence negligible. The threshold parameter β determines whether the systems is excitable (β > 1) or exhibits self-sustained perodic oscillations (β < 1). In the following we consider β in the excitable regime. Fig. 1 shows a schematic phase portrait of a spatially homogeneous system in the excitable regime (β > 1) with the cubic activator nullcline and the vertical inhibitor nullcline (solid red and green lines). The system has a single fixed point (u s , v s ), which is stable for β > 1 and located on the left branch of the cubic nullcline. At β = 1 a supercritical Hopf bifurcation of a limit cycle occurs, and the fixed point becomes unstable and is shifted to the middle branch of the nullcline for β < 1.
The excitable behavior of the system is crucially determined by the cubic nonlinearity of the activator equation and the separation of time-scales between the two variables: When the system is perturbed by a sufficiently large (super-threshold) external stimulus, which can be regarded as setting the initial condition, the system undergoes a large excursion in phase space (spiking). Starting from its initial condition, the system performs, due to the strong timescale separation ε 1, a fast transition to the stable right branch of the activator nullcline. After that, it travels slowly upwards approximately along this nullcline, until the phase points jumps back to the left branch, and returns, along the left branch of the nullcline, slowly downwards to the fixed point (recovery phase). Without further external stimulation the system remains in the stable fixed point (rest state).
The threshold-like behavior of the FitzHugh-Nagumo system is associated with the canard-like trajectory (dotted black line in Fig.( 1)), which is the trajectory passing through the local maximum of the cubic nullcline, and which is often referred to as the threshold of the FitzHugh-Nagumo system. The region around the canard-like trajectory is extremely sensitive to initial conditions: For initial conditions only slightly below the canard-like trajectory the systems will perform a large excursion in phase space, whereas for initial conditions only slightly above the canard-like trajectory the excursion will be small (sub-threshold excitation). In principle, the transition from small to large amplitude excitation is continuous; in fact, however, phase space excursions of intermediate amplitude are very rare. Correspondingly, small sub-threshold stimulations (dash-dotted black line in Fig.( 1)) will result in fast relaxation, while super-threshold stimulations (dashed black line in Fig.( 1)) induce a full excursion in phase space, corresponding to a characteristic spike in the time evolution of the u-variable.
In the following we consider spatially inhomogeneous solutions, i.e., waves or wave segments in two spatial dimensions (2D). Depending upon the set of parameters (ε, β) there exist different wave solutions [34,5]. Here we focus on localized wave segments which may either shrink or expand, as they propagate, or in the limit case, remain unchanged in size and shape, in which case they are called particle-like waves or dissipative solitons.
As we are interested in stationary propagating wave segments, Eqs. ( 1)-( 2) can be written as where ξ = x + ct is the co-moving coordinate with propagation velocity c.
In a co-moving frame, the localized critical structure is related to a saddle-point with a single unstable eigenvector (one-dimensional unstable manifold) in phase space. The curve representing this solution in a parameter plane of the bifurcation diagram is called the rotor boundary ∂R. In phase space, the stable manifold of states on ∂R separates the attractor of a spiral wave (spatially non-confined) and the stable, spatially uniform steady state. Thus, when the form of this dissipative soliton is perturbed, it either grows to a spiral wave, or shrinks and disappears. Perturbations above the critical size of the dissipative soliton grow, while perturbations below that critical nucleus shrink to the stable uniform state, i.e., the wave segments retracts. The internal cortical control of such dissipative solitons may be viewed as a strategy of the cortex to avoid re-entrant spiral waves, e.g., in migraine. Changing the nucleation size of this critical structure changes the susceptibility to pathological conditions such as spreading depression.
In Fig. 2, we show the rotor boundary ∂R (black dashed) in a schematic bifurcation diagram of wavesize S as a function of the threshold parameter β. It separates the weakly excitable parameter regime (perturbations grow to a spiral wave) from the subexcitable parameter regime (wave segments in 2D below the critical size shrink in length, while in spatially one-dimensional (1D) systems wave propagation is stable). There exists another boundary ∂R ∞ (dash-dotted vertical line), independent of size S, to the right of which all wave segments retract (corresponding to infinitely large critical size). Furthermore, the propagation boundary ∂P 1D is shown, which separates the subexcitable parameter regime from the nonexcitable regime. In the nonexcitable parameter regime, perturbations shrink also in width and wave segments collapse, i.e., even in spatially one-dimensional systems no wave propagation is possible.
At this point we would like to remark that the critical nucleus (dissipative soliton), which has the dynamic signature of a saddle-point, can be stabilized by an internal feedback control loop which controls the excitation threshold β in Eq. ( 2) that is, where K is the control strength, and the size of the wave segment S represents a measure of the active area occupied by the wave segment.
Eq. (5) defines a control line (red solid line with arrows) in the (β, S) phase diagram in Fig. 2. As the temporal evolution of a perturbation in a controlled system is confined to the control line, it asymptotically approaches a stable wave segment (β * , S * ) (if perturbed with convenient initial conditions). This follows from Fig. 2 since wave segments above ∂R, i.e., β < β * , grow in size, while wave segments below ∂R, i.e., β > β * , shrink.
The size of the wave segment S can be calculated via an integral over the excited area, e.g., If the stabilized solution (β * , S * ) is reached, S(t) becomes stationary. Thus, although the control is invasive, it does not produce new solution branches, and the stabilized critical nuclei of the controlled system are unstable solutions of the system without control at β = β(t −→ ∞). Hence, it is adequate to define the size of the wave segment as the area where the activator concentration u(r, t) is larger than zero with the Heaviside function Θ.
In the following simulations we will apply an internal feedback mechanism as in Eq. (5) with the wavesize S(t) from Eq. (7) to stabilize the critical nucleus. Furthermore, we will study the influence of the curvature of the excitable medium on the stability of localized waves. Hence the Laplace operator ∇ 2 must be replaced with the Laplace-Beltrami operator ∆ LB [35] for surfaces given in curvilinear coordinates α i , i = 1, 2,: where g = Det G and G with the matrix elements g ij is the metric tensor of the parametrization, see Appendix.
The surface of a torus in the Euclidian space R 3 can be described by the parametrization (θ, ϕ) of the position vectors The geometrical meaning of the major curvature radius R and the minor curvature radius r as well as the angles θ and ϕ is visualized in Fig. 3. The Laplace-Beltrami operator in torus coordinates reads We investigate sections of a torus with von Neumann boundary conditions (no flux boundary) on the equatorial section (at θ = 0 and θ = π) and periodic boundary conditions in the direction of the azimuthal angle ϕ. This restricts all traveling wave solutions as they have to obey the symmetries defined by these boundary conditions, i.e., the center of mass of the critical nucleus is pinned either on the outside or inside of the torus.
The FHN model with diffusion and internal feedback control Eq (5) including the wavesize S(t) calculated with Eq (7) ∂u is numerically solved using the explicit Euler method. With the discretization of space and time where x stands for θ or ϕ, respectively, the time derivative is calculated as with u n j := u(t n , x j ). The Laplace-Beltrami operator Eq. (8) contains derivatives of first and second order with respect to θ and ϕ. The derivative of first order is solved using a forward-backward-Euler-algorithm The derivatives of second order are solved using the Euler method, first backward followed by forward As initial condition, the activator concentration u is set equal to u s + 2 (which corresponds to a supra-threshold excitation in Fig. 1) in a rectangular area, and the inhibitor concentration v is set equal to v s + 1.5 in a rectangular area shifted relative to the activator stimulus, in order to determine the propagation direction of the wave segment. Outside the rectangle, the initial condition is u = u s , v = v s , with the activator and inhibitor fixed point values u s , v s . In order to numerically obtain critical nuclei of smaller or larger size, the size of the rectangular area is varied.

Overview of wave solutions on a torus
The main results are twofold. First, by investigating excitation waves on a torus, we show that the Gaussian curvature of the excitable medium changes the nucleation threshold in a systematic way. Second, and more surprisingly, we observe that a curved and major curvature radius R = 80 2π ; 2: solutions on a torus with minor curvature radius r = 20 2π and major curvature radius R = 40 2π ; stable ring wave solutions (red solid lines) with points of excitation block, i.e., propagation suppression (red asterisks); unstable inside critical nucleus (blue dash-dotted lines); unstable outside critical nucleus (green dashed lines); stable stationary and stable oscillating localized wave segment on the torus outside (green solid lines). Feedback Eq. (5) is applied to stabilize the states on the dashed and dash-dotted curves. medium can even induce a change of stability. Unstable critical nuclei are transformed into stable propagating localized wave segments.
We analyze the nucleation of excitation waves in reaction-diffusion media on curved 2D surfaces, specifically on tori. A torus has positive and negative Gaussian curvature on the outside (θ = 0) and the inside (θ = π), respectively, and a continuous transition in between with vanishing Gaussian curvature on the top (θ = π 2 ) and bottom (θ = 3π 2 ), see Fig. 3. In general, a torus has, in contrast to a sphere, not only locally varying and even negative Gaussian curvature, but a torus also admits a global isothermal coordinate system, called toroidal coordinates ((θ i ,φ), see Appendix), that is, coordinates where the metric is locally conformal to the Euclidean metric. Therefore an intuitive understanding of some of our results can be based on the particularly simple form of the Laplace- Beltrami operator, the "diffusion operator", in these coordinates.
On tori, a stable solution besides the spatially homogeneous steady state are ringshaped localized traveling wave solutions. Ring-shaped traveling waves have been analyzed, and in particular their critical properties have been discussed, namely, wave fronts with sufficiently large geodesic curvature break up on the torus inside [36]. We reconsider these travelling waves in order to compare them with the dynamics of critical nuclei. The stable manifold of a critical nucleus separates the attractor of a ring-shaped traveling wave and the spatially uniform steady state.
We restrict our study to nucleation of waves propagating strictly in the direction of the azimuthal angle ϕ (see Sect. ( 2)) and, furthermore, the center of mass of the nucleation is pinned either on the outside or inside of the torus, i.e., the locations where the extreme values of the Gaussian curvature occur. In the following, we will simply refer to these solutions as inside or outside critical nuclei or, if stabilized, inside or outside traveling wave segments. These solutions are the symmetric solutions with respect to the equatorial section of the torus. Note that there exist also asymmetric solutions on tori that we have not studied. Furthermore, note that these two symmetric critical nuclei on the outside and inside assume the shape of a localized wave segment where the open ends extend in the direction of θ (perpendicular to the propagation direction), that is, into regions of decreasing and increasing Gaussian curvature, respectively. Our results are mainly explained by this gradient in the Gaussian curvature and not by the absolute value of the Gaussian curvature.
The results are displayed in two bifurcation diagrams. First, the same bifurcation diagram as already introduced in Sect. 2 to define the regimes of excitability (weakly, sub-, and nonexcitable, see Fig. 2 in Sect. 2) is shown in Fig 4. The size S of the critical nucleus, see Eq. (7), is plotted versus the threshold parameter β of the local dynamics of the FHN system Eq. (1)- (2). The reference branch of the critical nucleus from simulation on a flat medium is now labeled "flat" in Fig 4 (black dashed). It separates the weakly excitable regime (to the left, decreasing β) from the subexcitable regime (to the right, increasing β), which ends at β = ∂P 1D , where the nonexcitable regime is reached.
In Fig 4, we show further solution branches simulated on two different tori. The torus labeled 1 has lower absolute values of Gaussian curvature than the torus labeled 2, since the latter torus has a two times smaller value of its major curvature radius R. For each torus, we have a branch of the ring-shaped traveling wave solution (red solid). Furthermore, for each torus, we have one branch of the inside (blue dashdotted) and outside (green dashed) critical nucleus. The states on the curves of the critical nucleus (dashed or dash-dotted) are stabilized by applying an appropriate global feedback Eq. (5) with suitably chosen β 0 and K such that the respective state (β * , S * ) is at the intersection with the line given by Eq. (5). In addition, on the torus outside, we find stable wave segments and stable oscillating waves (green solid), see Sect. 3.4 below.
Second, Fig. 5 is a bifurcation diagram, where the propagation velocity c, see Eq. (3)-(4), is plotted versus the threshold parameter β. The reference branches are on the one hand the propagation velocities of the stable fast and the unstable slow wave solution in spatially one-dimensional systems with diffusion in one spatial direction only, [37]), (grey dashed), and, on the other hand, a critical velocity c min (black solid), below which stable wave propagation cannot be obtained [36].
In For the outside stable wave segments, we show the propagation velocity c at the center of mass (green solid) and, after the bifurcation into two branches with decreasing threshold parameter β, the maximum and minimum propagation velocity of the stable oscillating wave segment. In addition, for the outside critical nucleus (with propagation velocity c o ), we plot a "hypothetical" branch (green dashed) that shows the propagation velocity c i = c o R−r R+r which a point of this wave segment would have on the torus inside if it existed there.

Ring wave break-up at saddle-node bifurcation
First we focus on the break-up of ring-shaped traveling waves on tori [36]. The ringshaped traveling wave solution is a stable solution of Eqs. (1)-(2) shown in Fig. 6. Thus the ring waves can be concieved as homoclinic solutions of the related ordinary differential equations (3),(4) in the co-moving frame. Ring waves have negative geodesic curvature on the torus inside and positive geodesic curvature on the torus outside, see Fig. 6(a). Thus, compared to 1D pulses (or infinitely extended wavefronts on a flat surface, respectively), ring waves propagate more slowly on the torus inside and faster on the torus outside, see Fig. 5 (red solid). If the propagation velocity falls below a critical value c min , the ring wave breaks up on the torus inside [36], see Fig. 6(b). This excitation block is marked by an asterisk in Figs. 4 and 5. Below the minimal velocity c min , stable wave propagation cannot be obtained. For 1D waves, it is known that the propagation failure is due to the coalescence of the homoclinic orbits of the fast and the slow wave (pulse) solution of the ODE problem Eq. (3)-(4) [37]. In the related PDE problem Eqs. (1)-(2), the propagation boundary is a saddle-node bifurcation point of the stable fast wave branch and the unstable slow wave branch. In Fig. 5, the fast wave branch and the slow wave branch of 1D traveling wave solutions are shown (upper and lower dashed grey line). At the propagation boundary ∂P 1D , they meet in a saddle-node bifurcation. Also in curved 2D media, the propagation failure is due to a saddle-node bifurcation, where the fast wave branch collides with the slow wave branch in a saddle-node bifurcation. In the 1D limit R −→ ∞, the threshold β, at which the ring waves break up, converges to the propagation boundary ∂P 1D . This is shown in Fig. 7, where the lines show the propagation failure in the (R, β) parameter space on two different tori; the upper line (dashed blue) is computed on a less curved torus with lower absolute values of Gaussian curvature compared to the lower line (dash-dotted green). The minimum velocity c min , below which stable wave propagation is not possible, can be calculated as (see Appendix) where ε cr is the critical time separation parameter, where the homoclinic orbits of the 1D fast and slow wave (pulse) solution of the ODE Eq. dynamics (ε and β) and the Gaussian curvature Γ of the torus. An increase of β or ε causes a deceleration of the ring wave on the torus inside. Also an increase of the Gaussian curvature Γ of a torus causes an increase of the absolute values of the geodesic curvature of the ring wave, which results in a deceleration of the ring wave on the torus inside. Thus, on more strongly curved tori (smaller R), the ring wave breaks up at smaller threshold β, see Fig. 7.

Curvature-induced changes of nucleation
Next, we analyze the nucleation of excitation waves on the torus inside and outside, respectively.
Torus inside The inside branches of the critical nucleus (blue dashed) in Fig. 4 are to the right of the reference curve (rotor boundary ∂R on flat surfaces). The larger the size S of the critical nucleus is, the stronger is the shift towards larger threshold β. On the more strongly curved torus (torus 2), the branch of the critical nucleus is shifted further. Thus, on the torus inside, there exist critical nuclei in the subexcitable parameter regime, where on flat surfaces all wave segments retract and vanish, cf. Fig. 2. A qualitiative explanation of this behaviour can be given by the relation of the Gaussian curvature Γ at the center of mass of the critical nucleus (θ = π) and at the open ends of the critical nucleus. Mathematically, the Gaussian curvature is described by the Laplace-Beltrami operator Eq. (8) in torus coordinates [35]. A torus admits a global isothermal orthogonal coordinate system, so-called toroidal coordinates (θ i ,φ), that is, coordinates, where the metric is locally conformal to the Euclidean metric. The Laplace-Beltrami operator Eq. (8) given in toroidal coordinates is where a = (R 2 − r 2 ) 1 2 is a measure for the scaling of the space, η = arcoth[R/(R 2 − r 2 ) 1 2 ] is a measure for the relation between the major radius R and the minor radius r and ϕ = ϕ sinh η. The derivation can be found in the Appendix. Introducing an effective coupling strength C(θ i ) = (cosh η − cos θ) 2 /a 2 , a torus can mathematically be interpreted as a flat medium with a spatial coupling being a function only of the location θ (θ i can be expressed in terms of θ, see Appendix), i.e.D(θ) = DC(θ). This is similar to effective geometric potentials arising from curved surfaces in hard condensed-matter physics [38]. The coupling strength C(θ) is strictly monotonically increasing from the torus outside (θ = 0) to the torus inside (θ = π), see Torus outside On the torus outside, we find that, under certain conditions, unstable critical nuclei bifurcate into stable propagating wave segments (green solid line in Fig. 4).  Furthermore, we find stable oscillating wave segments, whose size oscillates periodically in a self-sustained way. This striking bifurcation pattern will be explained in Sect. 3.4. The outside branches (green) in Fig. 4 are to the left of the flat reference branch ∂R (black dashed). On more strongly curved tori, the branch of the critical nucleus is further shifted. The coupling strength relation between the torus inside and outside also explains this behaviour. As the coupling strength C(θ) at the center of mass of the outside critical nucleus is smaller than at the open ends, the resultant diffusion perpendicular to the propagation direction is directed towards the center of mass, which enhances the retraction of the open ends.
On the torus outside, critical nuclei with increasing size S are found at decreasing threshold β. This is distinct from the inside nucleation branch and the flat reference branch: the larger the size S of the critical nucleus is, the larger is the difference between the coupling strength at the center of mass and the coupling strength at the open ends. Thus, larger critical nuclei at the torus outside are shifted to smaller threshold β, whereas on the torus inside larger critical nuclei are shifted to larger threshold β.
Critical nuclei with small size S extend over an area with almost constant coupling strength C, see Fig. 8. This supports the assumption that the branches of the inside and outside critical nucleus with small wavesize S in Fig. 4 lie close to the flat reference branch. This implies that the outside nucleation branch (green solid in Fig. 4) at small wavesize S terminates in a saddle-node bifurcation, and the unstable outside nucleation branch coalesces with the flat reference branch; this could, however, not be resolved numerically.

Curvature-induced stabilisation
Depending upon the excitation parameter β, different space-time patterns occur (Fig. 9). Localized wave segments may either grow to become stable ring waves ( Fig. 9(a)), or they may shrink and vanish ( Fig. 9(d)). Additionally, as an effect of the curved surface, on the torus outside, we find stable propagating localized wave segments, see Fig. 9(c). Furthermore, we find stable oscillating wave segments, whose size oscillates periodically, see Fig. 9(b). In Fig. 10, we show the activator profile of a stable wave segment propagating with a stationary profile, and in Figs. 11(a) and (b) we show snapshots of the activator profile of an oscillating wave segment at its minimum size S and at its maximum size S, respectively. The existence of stable wave segments on surfaces with positive Gaussian curvature can be qualitatively explained with the help of the space-dependent effective coupling strength C, as discussed in Sect. 3.3. The open ends of a stable wave segment on the torus outside lie in an area of the torus where the coupling strength C is larger than the coupling strength at θ = 0, where the center of mass of the wave segment is located. Thus the resultant effective diffusion perpendicular to the propagation direction caused by curvature is directed towards the center of mass of the wave segment. The larger the size S of the perturbation is, the stronger is this effect. At the same time, in the excitable parameter regime (see Fig. 2), small perturbations grow in length. If these two effects are balanced, we find stable propagating wave segments. In Fig. 4, we show the branch of stable wave segments (green solid). Perturbations with size S larger than the stable wave segments (and smaller than ring waves) shrink, as the difference in effective coupling strength between the center of mass and the open ends is large. Perturbations with size S smaller than the stable wave segments and larger than the small outside critical nucleus (which is not shown in Fig. 4 but supposed to lie close to the flat reference branch, see Sect. 3.3) grow, as the difference in coupling strength between the center of mass and the open ends is small.
In Fig. 5, we show the branch of the propagation velocity c at the center of mass of the stable wave segments and the stable oscillating wave segments (green solid). Furthermore, we show a hypothetical branch, the related propagation velocity at the torus inside (green dashed).
It is impossible that the stable wave segments grow to ring-shaped traveling waves (without enlarging their geodesic curvature), as the hypothetical propagation velocity at the torus inside is smaller than the minimal velocity c min (black solid line in Fig. 5).
For decreasing threshold β, the hypothetical propagation velocity at the torus inside of the stable outside wave solution (green dashed) accelerates, whereas the minimum velocity c min (black solid) slows down. At the intersection point of these two branches, the stable wave solution bifurcates into a stable oscillating wave segment and an unstable critical nucleus.

Conclusions
Spreading depression (SD) is a pathological dysfunction of brain activity that occurs, e.g., during migraine or stroke. It appears as spatially localized wave segments propagating in the cortex [5]. Despite substantial progress in the understanding of SD, the inherent processes are still incompletely known. Here, we study the influence of the curvature of the cortex on SD. For this purpose, nucleation and propagation of spatially localized reaction-diffusion waves are investigated on the surface of a torus. These unstable structures, which represent critical nuclei, are stabilized by an internal feedback control Eq. 5. We confine attention to the case that the center of mass of the critical nuclei is pinned on the torus inside or outside, respectively.
We have shown that negative Gaussian curvature (torus inside) causes a shift of the nucleation branch to larger threshold β, i.e., there exist critical nuclei in a parameter regime that is subexcitable on flat surfaces. In view of SD waves in the cortex, this might indicate that SD is more likely to be initiated in areas with negative Gaussian curvature.
In addition, we have made the surprising discovery that curvature can induce a change of stability, i.e., on the torus outside we have found stable propagating localized wave segments, as well as wave segments periodically oscillating in size.
These findings are qualitatively explained by an effective coupling strength, which can be found as an equivalent mathematical description of the Gaussian curvature on surfaces that admit isothermal coordinates.
Furthermore, we have reviewed the behaviour of ring-shaped wave solutions (autowaves), which were first described by V. A. Davydov in 2003 [36]. In the context of critical propagation effects, we have confirmed that the propagation boundary of ring-shaped waves, constituting break-up on the torus inside, is caused by a saddle-node bifurcation, where the fast wave branch collides with the slow wave branch. Thereby, the propagation velocity of ring-shaped waves is compared to a semianalytically calculated minimum velocity c min Eq. 17.

Minimum propagation velocity
It is well known that wave propagation becomes impossible if the propagation velocity of a traveling wave falls below a critical value c min [39]. Here, we consider the case of a traveling wave on a curved surface, which leads to slowing down of the wave by negative geodesic curvature of the wave front. The geodesic curvature of a wavefront is the curvature of the front projected in its tangential plane. This should not be confused with the Gaussian curvature Γ of a surface, which is the product of the two principal curvatures, i.e., Γ = cos θ r(R + r cos θ) (19) in common torus coordinates (θ, ϕ). The FHN model (1), (2) in one spatial dimension for 1 < β < √ 3 and sufficiently small ε, has a stable fast wave solution and an unstable slow wave solution corresponding to homoclinic orbits of the related ODE problem Eqs. (3),(4), see [37]. There exists a critical line in the (ε, β) space, at which the fast wave branch is connected to the slow wave branch. For values of β or ε, respectively, above this critical line, propagation of traveling waves is not possible. The values of ε on this critical line are denoted by ε cr , the corresponding critical propagation velocity is c cr .
Here we show the derivation of the analytical dependency between the minimum velocity c min and the critical time scale separation ε cr and the critical propagation velocity c cr .
As described in [39], the propagation of slightly curved fronts (curvature radius R rising front width L) can be approximated by where the FHN model has been transformed to the co-moving coordinate ξ = x + ct, with propagation velocity c. This can be written as Introducing rescaled parameters c * and ε * yields which has the same form as the FHN model in one spatial dimension. Thus, for c * = c cr and ε * = ε cr , the homoclinic solution of Eqs. (26), (27) corresponds to the connection between the fast wave branch and the slow wave branch. Hence, the minimal velocity c min of curved fronts can be derived from Eq. (25) by setting c * = c cr and ε * = ε cr . This yields c min = ε ε cr c cr .
Toroidal coordinates A parametrization f : {α i } → {x j } gives the Laplace-Beltrami operator in curvilinear coordinates where G is the metric tensor with matrix elements g ik , which is the product of the transposed Jacobian matrix of f multiplied with the Jacobian matrix of f , and g = Det G, see [35]. The single components of the metrical tensor thus are the scalar product A parametrization f is isothermal, if the derived coordinate system is orthogonal and conformal. In two spatial dimensions, a parametrization is orthogonal, if the scalar product of the basis vectors for i = k equals zero, ∂f ∂α i | ∂f ∂α k = 0.
The condition for conformal mapping is ∂f see [35]. This yields the following form of the Laplace-Beltrami operator To derive a global isothermal coordinate system for the surface of a torus, we start from the parametrization [40] (θ i , ϕ) →    a sinh η cos ϕ cosh η−cos θ i a sinh η sin ϕ cosh η−cos θ i a sin θ i cosh η−cos θ i where a > 0 is a scaling factor of space and η > 0 is a measure for the ratio of major curvature radius R and minor curvature radius r. As can easily be proved, these are orthogonal coordinates. As g θ i θ i = a 2 (cosh η − cos θ i ) 2 and g ϕϕ = a 2 sinh 2 η (cosh η − cos θ i ) 2 , g θ i θ i = g ϕϕ , thus the parametrization (35) is not conformal. Introducing the variablẽ Thus the Laplace-Beltrami operator in isothermal torus coordinates reads To obtain the dependencies of a and η upon the major curvature radius R and the minor curvature radius r, which are parameters of the common parametrization (θ, ϕ) →    (R + r cos θ) cos ϕ (R + r cos θ) sin ϕ r sin θ one needs to compare Eq. (38) with the isothermal parametrization a sinh η cos(φ sinh η ) cosh η−cos θ i a sinh η sin(φ sinh η ) cosh η−cos θ i a sin θ i cosh η−cos θ i A necessary and sufficient condition that a point from the domain of definition of the parametrization Eq. (38) lies on the twodimensional surface of a torus in the Euclidian R 3 is In toroidal coordinates Eq. (39) yields x 2 + y 2 + z 2 − 2a cosh η sinh η x 2 + y 2 + a 2 = 0.
To derive the dependency between θ i and θ, the expressions x 2 + y 2 − R of both coordinate systems are compared. This yields r cos θ = a sinh η cosh η − cos θ i − R.