Excitable media store and transfer complicated information via topological defect motion

Excitable media are prevalent models for describing interesting effects in physical, chemical, and biological systems such as pattern formation, chaos, and wave propagation. In this manuscript, we propose a spatially extended variant of the FitzHugh–Nagumo model that exhibits new effects. In this excitable medium, waves of new kinds propagate. We show that the time evolution of the medium state at the wavefronts is determined by complicated attractors which can be chaotic. The dimension of these attractors can be large and we can control the attractor structure by initial data and a few parameters. These waves are capable transfer complicated information given by a Turing machine or associative memory. We show that these waves are capable to perform cell differentiation creating complicated patterns.


Introduction
Excitable media are popular models for the exploration of physical, chemical, and biological systems. In this paper, we propose a modification of the FitzHugh-Nagumo (FN) model that exhibits new intriguing phenomena. This model describes an excitable media where the waves of a new kind propagate. We show that the time evolution of the medium state at the wavefronts is determined by complicated attractors that can be chaotic. The dimension of these attractors can be made large and increased to infinity depending on the model parameters. We can control the wavefront dynamics structure by initial data and a few parameters. These waves are capable transfer complicated information given by a Turing machine or associative memory. The physical mechanism of the formation of the waves is based on an interaction between topological defects (kinks). Recently we have had an understanding that topological defect motions in an active medium are important for morphogenesis [1][2][3][4], and we show that these waves are capable to perform cell differentiation creating complicated patterns.
First, we outline our model which we consider as a spatially extended variant of the FN model. Similar to the FN model, the model involves cubic nonlinearities, but we extend it by including small spatial gradients. Then the first equation of the model can be considered as a perturbed scalar Ginzburg-Landau (GL) equation for a scalar order parameter u with a small gradient term. That equation describes bistability and spontaneous layered patterning. The GL equation simulates a trigger mechanism that in real biological systems is generated by positive feedback loops in gene regulation networks [5]. Moreover, we implement diffusion effects in the second equation of the FN model and we complement the two equations for components u, v by the linear wave equation for a third variable z that would be interpreted as mechanical deformation. The equation for z is coupled with two first ones for variables u, v via a quadratic term. Now we outline the physical mechanism of wave generation. Let ϵ 2 /2 be the coefficient at the gradient term in the corresponding GL energy and let ϵ > 0 be small. It is well known that such a singularly perturbed GL equation has asymptotical solutions describing kink chains where ith kink is localized at some positions slowly evolving in time [6]. Kinks are narrow topological defects (of the width O(ϵ)) with the charges (−1) i describing a symmetry breaking and a layered pattern formation: a separation of the entire domain on subdomains along x-axis, where u ≈ ±1. Note that the direct interaction between kinks is exponentially small (it has the order O(exp(−c 1 /ϵ))). Therefore such asymptotical solutions are correct within an exponentially long time O(exp(−c 2 /ϵ)) while kinks are well-separated [6]. We consider a simple perturbation of the scalar Ginzburg-Landau equation by the second variable v. The time evolution of v is defined by the reaction-diffusion dynamics where the order parameter u is involved via the quadratic term zu x . The time evolution of z is governed by a linear wave equation. We suppose that the v-component diffuses fast. The u-kinks interact with the fast reagent v and in turn that reagent v acts on the u-kinks. We thus have a feedback loop that leads to a small non-local kink interaction via the intermediate v-component. This interaction is small nonetheless it is much stronger than the direct local interaction between kinks that is exponentially small as ϵ → 0. Note that our system for u and v reminds the model studied in [7] where the formation of layered patterns is studied.
Actually, our mechanism of the kink motion and layered pattern formation is like to [7], however, there is a key difference: our system for u, v is open, and, due to the presence of the term zu x there exist no Lyapunov functions decreasing along trajectories and making the sense of energy. This property is very important. In fact, we show that under an appropriate choice of system parameters the dynamics of the kink coordinates X i can be described by the Hopfield system of ordinary differential equations with continuous-time and, in general, non-symmetric interactions. This non-symmetry is a consequence of the fact that our system is open. It is well known that such Hopfield systems exhibit a remarkable universality property [8]: they can generate any structurally stable (hyperbolic) dynamics. Such dynamics may be chaotic (the best-known examples are given by Anosov flows and Smale horseshoes [9,10]). Following [11,12] we can use this chaos to simulate Turing machines and we apply it to cell differentiation using the results from [13].
Thus, combining reaction, diffusion, and wave effects we obtain an excitable media that can generate waves that can be considered as moving neural nets consisting of interacting narrow fronts (kinks). The evolution of the coordinates that define the localization of those fronts is governed by a dynamical system. The key point is that we can control the attractors of these dynamical systems by positional information stored in spatially distributed initial data for z and by the choice of a few parameters. These attractors may be chaotic and of high dimension. Information on the attractor structure is contained in initial data for z-component. So, our excitable media transform positional information stored in space into dynamic information and complex temporal behavior. We think that such waves can transfer adaptive information in cell colonies and multicellular organisms.
Using results [13] we will show that by such waves one obtain any layered patterns of cell differentiation (see Section 7). The idea that topological defects can serve as ''topological morphogens'' has received attention last decade and studied theoretically and experimentally (see, for example, [1][2][3][4]).
The paper is organized as follows. In the next section, we describe the model. In Section 3, we state the property of Universal Dynamical Approximation and by this mathematical formalism, we state the main results in Section 4. The subsequent sections contain an outline of the proof, construction of asymptotical solutions, application to morphogenesis, and concluding remarks.

Model
The model consists of a scalar perturbed Ginzburg-Landau equation for an order parameter u, a reaction-diffusion equation, and a hyperbolic equation for z: Here u = u(x, y, t) and v(x, y, t) are unknown functions defined on Ω × {t ≥ 0}, Ω is the strip (−∞, ∞) × [0, 1] ⊂ R 2 , κ > 0 and ϵ > 0 are small parameters. To simplify the problem and provide existence of solutions for all times t > 0 we set the 2π -periodic boundary conditions for v, u and z supposing that initial data satisfy these periodicity conditions.
At the horizontal boundaries y = 0 and y = 1 we set the Dirichlet conditions for v: and the zero Neumann condition for u u y (x, y, t) The initial conditions are given by smooth functions u 0 , v 0 and z 0 , for example, and similarly for u, v. The function z 0 plays a key role in the control of wave large time dynamics.
This system reminds the celebrated FitzHugh-Nagumo (FN) model [14,15] but it is extended by spatial gradients and the additional equation for a variable z that can describe mechanical effects critically important in biological active media [1,2]. So, this model takes into account basic mechanical, chemical, and physical effects. Note that our model is twodimensional which is important for the control of large time dynamics. To obtain analogous results in a one-dimensional case, we have to use a number of reagents replacing a single Eq. (1) by a reaction-diffusion system [13]. Note moreover that the first equation arises in the generalized Landau-de Gennes model applied for morphogenesis problems in [1] (if we remove the terms connected with cell polarities). To conclude this section, let us make some remarks. First, the model proposed can be simplified (although then the mathematical analysis is slightly more sophisticated). For example, one can remove the term κu x in Eq. (1).
The second remark concerns the generation of complicated dynamics without waves. The complicated attractors can be obtained in a shorted model consisting of Eqs. (1) and (2) with κ = 0. Then the coefficient z = z(x, y) can be considered as a parameter.

Universal dynamical approximation
To formulate the main results let us outline first the concept of Universal Dynamical Approximation (UDA) (that term is inspired by the work [16]). To explain this concept, remind that many neural networks such as multilayered perceptrons (MLP) enjoy the property of universal approximation: for each prescribed sufficiently smooth output function f defined on a compact domain D ⊂ R n of a Euclidean space R n and each ϵ > 0 we can adjust parameters of MLP in such a way that the output of f net (q) of the network is ϵ-close to f (q) for all entries q from the domain D [17].
The UDA concept generalizes universal approximation for the networks. Many evolution equations under reasonable boundary conditions define global semiflows in the appropriate functional phase spaces (see [18,19]). We denote such semiflows by S t P and they depend on the parameters P involved in equations, and boundary conditions. More formally, let us consider an evolution equation in a Banach space B depending on parameter P: Assume that for some P that equation generates a global semiflow S t . We obtain then a family F of global semiflows S t P where each semiflow depends on the parameter P.
Suppose for an integer n > 0 there is an appropriate value P n of the parameter P such that the corresponding global semiflow S t Pn has an n-dimensional normally hyperbolic locally invariant manifold M n embedded in our phase space B by a C 1 -smooth map defined on a ball B n ⊂ R n .

The restriction of semiflow S t
Pn to M n is defined then by a vector field Q on M n . Then we say that the family S t P realizes the vector field Q (that terminology is coined by P. Poláčik, see [20,21]).
The family S t P enjoys UDA if for each dimension n that family realizes a dense (in the norm C 1 (B n )) set of vector fields Q on the ball B n . Corollary 1. If the family S t P has UDA property, the Theorem on Persistence of Hyperbolic Sets implies that some semiflows S t P exhibit a chaotic large time behavior.
In other words, one can say that the UDA semiflows can simulate, by parameter variations, any finite-dimensional dynamics defined by a system of ordinary differential equations on the domain D ⊂ R n within any prescribed accuracy ϵ (in C 1 (D)-norm). This property implies that semiflows S t

Main results
Concluding the ideas presented above we formulate the following statements.
Theorem I. On dynamical complexity: under a choice of the parameters κ, ϵ, and initial data z 0 , u 0 (x) the kink dynamics of our model is defined by a time-continuous Hopfield system: whereX i are deviations of kink positions from some equilibrium valuesX i . This system has the property of universal dynamical approximation, where parameters are N, λ, h i and the entries K ij .
We can obtain any prescribed N, K ij , h i by a variation of ϵ > 0 and initial data for u and z.
This means that when we vary the model parameters, initial data and the kink number, kink coordinate dynamics can generate all possible kinds of structurally stable large-time behavior (up to topological equivalency). Since hyperbolic dynamics is persistent [10], kink motions generate all hyperbolic dynamics that may be chaotic [9,10].
The second result concerns with the shorted system (1) and (2). (1) and (2) under boundary conditions (5), (6) generates a family of global semiflows S t P , which has the UDA property, where parameters are ϵ > 0 and function z(x, y).
Note that averagesv(x, t) of v-pattern along y defined asv(x, t) = ∫ 1 0 v(x, y, t)dy can be described in a simple way as follows (up to small corrections vanishing as γ , ρ, ϵ → 0). These averages are sums of exponents So, if we know the velocities X i /dt of all kinks, we can restorev(x, t) for all x. To this end, we find the coefficients B m,j (t) using Eq. (11).
These results have interesting biological consequences. In Section 7, we show that positional information stored in a localized embryo domain can be transported in another domain and transformed into a developmental program.

Outline of the proof
The proof proceeds a few of steps. One can exclude z solving the wave equation and obtain a system for u, v with spatially heterogeneous coefficients. Analogous systems are well studied at physical level [7] and moreover there exists rigorous methods to analyze them [23]. We find kink chain solutions to the GL equation following [6]. These solutions are defined via kink coordinate X 1 , X 2 , . . . , X N . Under an appropriate parameter choice, x i slowly evolve in time. Although our system is two-dimensional, we can conserve planar kink structure by an appropriate choice of small parameters, and then kink dynamics can be completely defined by X i (t).
The next step is to derive dynamical equations for X i . It can be done quite rigorously in a standard way, and these equations capture all essential dynamics in our model while u remains close to kink chains. The equations for X i (t) form a system of coupled oscillators of a complicated form. Further, we show that under a special choice of z 0 that this system for X i can be simplified and reduced to a time-continuous Hopfield system (9), which describes small kink oscillations at certain positionsX i . This property of oscillation localization provides that the kinks remain well separated by distance δ ≫ ϵ for all times t. In the Hopfield system, deviationsX i of kink coordinates fromX i correspond to neuron states in neural network models. The parameters of this Hopfield system are the neuron (kink) number, the matrix K of the neuron interactions with the entries K ij and thresholds h i .
Next, there are two key points. The first one that in general K is non-symmetric thus there are no Lyapunov functions for the system (9). The second one is that we are capable to prove that we can control K: one can obtain any N × N square matrices K by a choice of initial data z 0 for z. This fact looks natural: to obtain a given square matrix K of size N × N, we should satisfy N 2 conditions by a countable set of unknown Fourier coefficients for z 0 (nonetheless, this fact needs a proof!). Just the second point on the control of K is verified, Theorems I and II at once follow from known results (for example, [22]).
We assume i.e., the minimal distance between the kinks is not too small: that minimal distance is much more than the characteristic diffusion length ϵ. We need condition (12) to obtain an asymptotical solution in the form of a kink chain.
The kink chains can be obtained by 2π -periodical in x functionsŪ N (x, X ). Inside narrow intervals I i,ϵ = ( , where s i = (−1) i are topological charges and i = 1, . . . , N. Outside of intervals I i,ϵ the functionŪ N (x, X ) is exponentially close to ±1 and this function is a smooth function of x. We can represent this kink chain as a sum of contributions associated with kinks: where U N,j is exponentially close to s j tanh(x − X j )/ϵ. For κ = 0 and v ≡ 0 we have a set of solutions u of (1) that have the form u =Ū N (x, X (t)) +ũ, whereũ is a small correction and the kink coordinates X j (t) evolve in time exponentially slowly: ) [6]. These solutions are correct while dist(X ) ≫ ϵ. The time evolution of X is a result of exponentially weak direct kink interaction. For κ ̸ = 0 we obtain the kink chain traveling with a constant speed κ: In the coming subsection, we consider the case v ̸ = 0, where there occurs a non-direct and non-local kink interaction via coupling with v-reagent.
The main idea in choice of γ is to conserve the planar structure of the kink fronts (otherwise it is impossible to describe kink chains by coordinates X i , and it is necessary to take into account the front curvature). The parameter κ > 0 should be small as well in order to obtain quasistationary solutions of Eq. (14). The condition c 0 exp(−c 1 ϵ −1/2 ) ≪ γ is necessary to ensure domination of non-local kink interaction via v-reagent with respect to direct kink one.
The subsequent statement follows works [6,8,23,26] with small modifications. Our first goal is to derive equations for X i (t). Eqs. for X i (t) can be obtained by a standard perturbation approach for small γ > 0 (see, for example, [6,26]). For sufficiently small γ > 0 one has u(x, t) = U N (x, X (t)) +ũ(x, t, γ ), (16) where U N (x, X ) is the kink chain (described above) andũ is a correction (the deformation of the kink chain form). The time evolution of X is governed by the equation 5 and small correctionsG i are uniformly bounded where s ∈ (0, 1], c, s are uniform in γ as γ → 0. To explain Eqs. (17) and (18), let us remind the construction from [23,26] that is the well known Lyapunov-Schmidt factorization. Let us introduce the standard convenient notation To obtain the dynamical equations for X i , we impose the condition ⟨ũ(·, ·, t), Ψ j (· − X i (t))⟩ = 0 (19) for each t and i. These equations define X uniquely for small ∥ũ∥, γ > 0 and bounded ∥v∥. Physically the conditions (19) mean that we split kink perturbations into kink shifts defined by X i and kink deformations orthogonal to the kink shifts.
For the correctionũ we then obtaiñ where L ϵ is the linear operator defined by and The spectrum of the operator L ϵ is well-studied [6]. Note that L ϵ is a self-adjoint operator of Schrödinger type that has a kernel consisting of N eigenfunctions, which are close to linear combinations of the kink Goldstone modes. We need the following Lemma (see [23]).
Lemma I. For X such that dist(X ) > δ 1 > 0 one has and if ⟨ψ, Ψ j ⟩ = 0 for all j then where all constants are uniform in ϵ as ϵ → 0.
For a proof see [23]. So, the spectrum of L ϵ consists of N exponentially small eigenvalues and all the remaining spectrum of L ϵ lies in the interval (−∞, −δ 0 ϵ 2 ), where δ 0 > 0 does not depend on ϵ > 0. This key spectral property implies that the system dynamics is defined by Goldstone modes via X i while the kink deformations remain small. This property also provides the stability of the kink solutions on exponentially large intervals I ϵ and allows us to solve Eq. (20). It can be done by the standard perturbative methods because for small γ , κ and ∥ũ∥ our model is weakly nonlinear while the linear part is stable due to condition (19).

Quasistationary solutions of (2)
Let us turn now to Eq. (2) for v. This equation is linear with respect to both v and u, u is a sum U N +ũ, whereũ is small. The function U N depends on time via slow variable X (t) and does not depend on t explicitly. Therefore, we can solve Eq. (2) by a simple idea: we can freeze X in the right-hand side of Eq. (2) assuming that X is just a parameter. The main contribution to v is given then by a function V N satisfying the equation We obtain 6 To obtain an asymptotic solution it is sufficient to note that for small ϵ > 0 the function dU N,j /dx = ϵ −1/2 Ψ j + O(exp(−cϵ −1 )) is a good approximation of δ-function (up to a constant uniform in ϵ > 0). Moreover, it is clear then that in Eqs. (25) z 0 (x, y) can be replaced by z 0 (X j , y). Let us denote by Γ m (x) the Green function of the one-dimensional boundary value problem satisfying the equation and the 2π -periodical boundary conditions in x.
Then we resolve (25) by the Fourier method that gives where we remove terms of the order O(ϵ) andẑ m (X j ) are the Fourier coefficients of z 0 (X j , y): where m are positive integers, b m = (1 − (−1) m )/(π m) (see Fig. 1).

Hopfield system
Using relations (24), Eqs. (17), (18), and removing small terms, we obtain evolution equations for kink coordinates: where Our goal is to reduce Eqs. (28) to the Hopfield system (9). It can be done by a special choice of z 0 (x, y), or, that is equivalent, ofẑ m (X j ). First we substitute formula (26) into (28). Then we have where The main idea to simplify the formula (29) for G i is as follows. Suppose that the kinks oscillate at certain fixed pointsX j , i.e., whereX i are new unknowns and ρ is a small parameter, which is independent of ϵ, γ , κ and defines the oscillation magnitude. Suppose temporarily that |X j | stay bounded as ρ → 0 (this assumption will be justified later). We can achieve such behavior of X -solutions of (28) under a special choice ofẑ 2m+1 (X j ). Positions of pointsX j may be arbitrary but the condition (12) must be satisfied. Namely, for x from a neighborhood ofX j we set where σ is a smooth sigmoidal function, for example, the Fermi function σ (z) = (1 + exp(−z)) −1 , and where M jm , S jm are unknown coefficients, which must be matched appropriately. Note that instead of the Fermi function we can take many other functions, for example, σ = sin.
We obtain then Fig. 1. The active medium described in the paper can generate waves, which can transfer a complicated time behavior. A cell colony that has created such a medium (or is immersed in it) can have an important selective advantage, watch the video [27]. For example, suppose that a cell, a member of that colony, finds a complex dynamical adaptive answer to an ecological challenge (the top row, the first cell). Then this answer can be transferred to other cells by the waves, and thus the whole colony obtains the ability to survive. Moreover, it is shown that the wavefront dynamics are defined by the Hopfield networks, so, those waves also may transport associative memory.
whereR ijm (X ) = O(ρ), are small corrections and Further, we use the following lemma.
Lemma II. For each N × N matrix K with entries K ij , i = 1, . . . ,N there exist a number M ≥ N and coefficients b jm such that Proof. For the unknown b jm we have a system of linear algebraic equations. For large m we have asymptotics for Γ m . Hence for sufficiently large M the matrix of our linear algebraic system contains a non-degenerate Vandermonde matrix as a submatrix thus that linear algebraic system is resolvable. 8 Using this lemma, we can choose S jm and M jm such that G i take the form where λ > 0 is a coefficient.

Control of dynamics for Hopfield system
Using (38) we obtain the system (9) for variablesX i .
It is easy to show that system (9) has a compact attractor. In fact, 0 ≤ σ (z) ≤ 1 thus that system implies the inequalities where |K| = max i,j |K ij |. These differential inequalities lead to the estimate The last estimate shows that system (9) has an absorbing set A = {X : |X i | < N|K|λ −1 }, thus it is dissipative and has a compact attractor. This result justifies our hypothesis on smallness of kink oscillationsX i at pointsX i and the transformation of general system (27) to the Hopfield system (9).
The following claim is proved in [8].
Lemma III. Dynamics defined by system (9) generates all finite dimensional hyperbolic dynamics (up to orbital topological equivalency) by variations of parameters K, N, λ and h.
In coming section, following [13], we show how this result can be applied for cell differentiation.

Applications to morphogenesis
In this section, we follow [13] however we propose here new ideas based on our model because this model essentially extends our possibilities with respect to model from [13] due to mechanical effects.
Multicellular organisms have evolved a diversity of cell types, typically, animal multicellular organisms contain tissues consisting of 100-150 different cell types. Our aim is to explain how waves can be used to describe cell differentiation and development of the whole embryo from a localized domain.
Let us note that usual traveling waves can be used to describe somitogenesis [28][29][30]. Somitogenesis is the process of somite formation which are bilaterally paired blocks that form along the anterior-posterior axis of the developing embryo in segmented animals. In vertebrates, somites give rise to different organs. The clock and wavefront (CWF) model describes the somite formation as a result of the oscillating expression of genes. The CWF model was proposed by Cooke and Zeeman [31] and developed by [30,32].
Let us outline the main ideas from [13]. Consider a layered one-dimensional (1D) pattern consisting of cells of two types, say, r and b. We decompose our ''organism'' into small domains of equal length occupied by cells. Then the pattern of cell differentiation can be considered as a string, for example, non-periodic one (rbrrrbbbr). Suppose that the cell type is determined by a morphogen operator M(u(·)) transforming patterns u(x) into patterns m(x 1 ), m(x 2 ), . . . m(x n ), where x j are cell centers and m(x j ) ∈ {r, b}.
To obtain a non-periodic pattern we need non-monotone concentration profiles that is not easy to obtain. Moreover, we would like to generate any prescribed patterns by a universal short model. To show how the system (1)-(3) really extends our possibility let us introduce a mathematical formalization of a complicated cell differentiation process. Let us consider strings s = (s 1 , s 2 , . . . , s n ) consisting of symbols s j ∈ {r, b} (extension to large alphabets is straightforward). The set of strings S D = {s (1) , s (2) , . . . , s (m) , the time moments t j = j∆T and pointsx 1 <x 2 ... <x n , where j = 1, . . . , n can be called a developmental program if the morphogen values take prescribed values at given time moments in given points: where j = 1, . . . , m. This formal definition is illustrated by Fig. 2.
Assertion DP. The shorted system (13)- (14) is capable to realize any developmental program by an appropriate morphogenetic operator M. This program can be encoded by initial kink positions X j and z 0 (X j , y).
This assertion is proved in SM.

Conclusions and discussions
In this paper, the active media are described where the generation of complex waves is possible. The waves are capable to transfer complicated (given by an attractor) information. This effect can be useful to explain classical morphogenesis effects such as the existence of embryo organizing centers. A developmental program created in a small domain can be transferred and released in another domain. Our model explains how it can be done by fundamental physical, chemical, and mechanical mechanisms. Namely, we use: (i) phase transitions via the Ginzburg-Landau model and kink formation, (ii) non-local kink interactions via an intermediate fast diffusing reagent, and (iii) linear mechanical waves. Note that kinks are the simplest topological defects, and in biological development, they can model segmentation, one of the most basic effects in morphogenesis.
It is clear that real biological active media are much more complicated (in particular, it can be described by the Landau-de Gennes model) so our extended FN system can be considered as a toy (but conceptual) model. An idea that topological defects (such as disclinations) may play an important role in morphogenesis and can be considered as ''topological morphogens'', became popular recently [1,2]. A mechanism consisting of elements (i)-(iii) permits to create of a complicated dynamical information stored in defect motion and look universal although of course for real media to proceed with rigorous mathematical analysis is a hard problem.
Note that a few genes is sufficient to encode the mechanism consisting of (i), (ii), and (iii). Therefore, it is natural to expect that such media could appear as a result of biological evolution. One can imagine, for example, such a situation. Consider a cell colony that must adapt to a new environment and develops products necessary for survival. It is clear that a colony, where it is possible to transfer complex adaptive innovations from one cell to another, has a clear selective advantage. This transmission can be done by means of the waves, opened in this paper, and these waves can not only transmit simple information, but they can also transfer complex behavior (which can be described by an attractor or a Turing machine, or a neural network with associative memory), similarly to human society.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Ivan Sudakow reports financial support was provided by National Science Foundation. Sergey Vakulenko reports financial support was provided by National Institutes of Health.