<figure><img src="image_1.png" /></figure>

Short communication

Excitable media store and transfer complicated information
via topological defect motion

<figure><img src="image_2.png" /></figure>

Check for
updates

Ivan Sudakow a,b,*, Sergey A. Vakulenko c,e, Dima Grigoriev d

ARTICLE INFO

Article history:

Received 17 February 2022

Received in revised form 8 July 2022

Accepted 17 July 2022

Available online 29 August 2022

Keywords:

Excitable media

Chaos

Waves

Cells

A B S T R A C T

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.

©2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ([http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)).

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

<footnote>ªSchool of Mathematics and Statistics, The Open University, Milton Keynes, MK7 6AA, UK</footnote>

<footnote>ᵇDepartment of Physics, University of Dayton, Dayton, OH 45469-2314, USA</footnote>

<footnote>ᶜInstitute of Problems in Mechanical Engineering, Russian Academy of Sciences, St. Petersburg 199178, Russia</footnote>

<footnote>ᵈCNRS, Mathématiques, Université de Lille, Villeneuve d'Ascq, 59655, France</footnote>

<footnote>ᵉSaint Petersburg Electrotechnical University LETI, St. Petersburg, 197022, Russia</footnote>

<footnote>* Corresponding author at: School of Mathematics and Statistics, The Open University, Milton Keynes, MK7 6AA, UK.
E-mail address: ivan.sudakow@open.ac.uk (I. Sudakow).</footnote>
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 $\epsilon^2/2$ be the coefficient at the gradient term in the corresponding GL energy and let $\epsilon > 0$ be small. It is well known that such a singularly perturbed GL equation has asymptotical solutions describing kink chains where $i$th kink is localized at some positions slowly evolving in time [6]. Kinks are narrow topological defects (of the width $O(\epsilon)$) 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 \approx \pm 1$. Note that the direct interaction between kinks is exponentially small (it has the order $O(\exp(-c_1/\epsilon))$). Therefore such asymptotical solutions are correct within an exponentially long time $O(\exp(-c_2/\epsilon))$ 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 $\epsilon \to 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-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.

## 2. 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$:

$$u_t = \frac{\epsilon^2}{2} \Delta u + u - u^3 - \kappa u_x + v, \quad (1)$$

$$v_t = \Delta v + z u_x, \quad (2)$$

$$z_{tt} - \kappa^2 \Delta z = 0. \quad (3)$$

Here $u = u(x, y, t)$ and $v(x, y, t)$ are unknown functions defined on $\Omega \times \{t \ge 0\}$, $\Omega$ is the strip $(-\infty, \infty) \times [0, 1] \subset \mathbb{R}^2$, $\kappa > 0$ and $\epsilon > 0$ are small parameters. To simplify the problem and provide existence of solutions for all times $t > 0$ we set the $2\pi$-periodic boundary conditions for $v, u$ and $z$

$$u, v, z(x, y, t) = u, v, z(x + 2\pi, y, t), \quad (4)$$

supposing that initial data satisfy these periodicity conditions.
At the horizontal boundaries $y = 0$ and $y = 1$ we set the Dirichlet conditions for $v$:

$$v(x, h, t) = v(x, 0, t) = 0, \qquad (5)$$

and the zero Neumann condition for $u$

$$u_y(x, y, t)|_{y=0,1} = 0. \qquad (6)$$

The initial conditions are given by smooth functions $u_0$, $v_0$ and $z_0$, for example,

$$z(x, y, 0) = z_0(x, y), \qquad (7)$$

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 two-dimensional 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 $\kappa 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 $\kappa = 0$. Then the coefficient $z = z(x, y)$ can be considered as a parameter.

### 3. 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 \subset \mathbb{R}^n$ of a Euclidean space $\mathbb{R}^n$ and each $\epsilon > 0$ we can adjust parameters of MLP in such a way that the output of $f_{\text{net}}(q)$ of the network is $\epsilon$-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_P^t$ 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$:

$$u_t = Au + F(u, P). \qquad (8)$$

Assume that for some $P$ that equation generates a global semiflow $S^t$. We obtain then a family $\mathcal{F}$ of global semiflows $S_P^t$ 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_{P_n}^t$ 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 \subset \mathbb{R}^n$.

The restriction of semiflow $S_P^t$ to $M_n$ is defined then by a vector field $Q$ on $M^n$. Then we say that the family $S_P^t$ realizes the vector field $Q$ (that terminology is coined by P. Poláčik, see [20,21]).

The family $S_P^t$ enjoys UDA if for each dimension $n$ that family realizes a dense (in the norm $C^1(\mathbb{B}^n)$) set of vector fields $Q$ on the ball $\mathbb{B}^n$.

**Corollary 1.** If the family $S_P^t$ has UDA property, the Theorem on Persistence of Hyperbolic Sets implies that some semiflows $S_P^t$ 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 \subset \mathbb{R}^n$ within any prescribed accuracy $\epsilon$ (in $C^1(D)$-norm). This property implies that semiflows $S_P^t$ can generate all structurally stable dynamics (up to orbital topological equivalency). Among the systems enjoying UDA, there are a number of fundamental ones: quasilinear parabolic equations [20,21], time-continuous and time recurrent neural networks [22], a large class of reaction-diffusion systems with heterogeneous sources [23], generalized Lotka-Volterra system [24], Oberbeck-Boussinesq model [25]. Also, the Euler equations on multidimensional manifolds exhibit similar properties [16]. Note that for time continuous and time recurrent neural networks and generalized Lotka-Volterra system the UDA property follows from Universal Approximation Theorem for MLP [22].
## 4. 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:

$$ \frac{d\tilde{X}_i}{dt} = \sum_{j=1}^{N} K_{ij}\sigma(\tilde{X}_j - h_j) - \lambda\tilde{X}_i, \quad (9) $$

where $\tilde{X}_i$ are deviations of kink positions from some equilibrium values $\tilde{X}_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).

**Theorem II.** System (1) and (2) under boundary conditions (5), (6) generates a family of global semiflows $S_p^f$, which has the UDA property, where parameters are $ε > 0$ and function $z(x, y)$.

Note that averages $\bar{v}(x, t)$ of $v$-pattern along $y$ defined as $\bar{v}(x, t) = \int_0^1 v(x, y, t)dy$ can be described in a simple way as follows (up to small corrections vanishing as $γ, ρ, ε \to 0$). These averages are sums of exponents

$$ \bar{v}(x, t) = \sum_{j=1}^{N} B_{m,j}(t) \exp(-m|x - \bar{X}_j|) \quad (10) $$

where $B_{m,j}(t)$ are coefficients such that

$$ \bar{v}(\bar{X}_i(t), t) = \text{const} \, dX_i/dt. \quad (11) $$

So, if we know the velocities $X_i/dt$ of all kinks, we can restore $\bar{v}(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.

## 5. 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, \dots, 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 positions $\bar{X}_i$. This property of oscillation localization provides that the kinks remain well separated by distance $δ \gg ε$ for all times $t$. In the Hopfield system, deviations $\tilde{X}_i$ of kink coordinates from $\bar{X}_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 \times 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 \times 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]).
# 6. Asymptotic solutions to the system (1)-(3) describing complex waves

## 6.1. Kink chains

Let us describe kink chains following [6]. Let first $\kappa = 0$. Let $X_j, j = 1, \dots, N$ be the coordinates of the kinks in the interval $(0, 2\pi)$. We suppose that $0 < X_1 < X_2 < \dots < X_N < 2\pi$. Let $dist(X) = \min_i X_{i+1} - X_i$, where formally $X_{N+1} = X_1$. We assume

$$ \epsilon \ll \delta < dist(X), \quad (12) $$

i.e., the minimal distance between the kinks is not too small: that minimal distance is much more than the characteristic diffusion length $\epsilon$. We need condition (12) to obtain an asymptotical solution in the form of a kink chain.

The kink chains can be obtained by $2\pi$-periodical in $x$ functions $\bar{U}_N(x, X)$. Inside narrow intervals $I_{i,\epsilon} = (X_i - \epsilon^{1/2}, X_i + \epsilon^{1/2})$, the functions $\bar{U}_N(x, X)$ have the form $s_i \tanh(\epsilon^{-1}(x - X_i))$, where $s_i = (-1)^i$ are topological charges and $i = 1, \dots, N$. Outside of intervals $I_{i,\epsilon}$ the function $\bar{U}_N(x, X)$ is exponentially close to $\pm 1$ and this function is a smooth function of $x$. We can represent this kink chain as a sum of contributions associated with kinks:

$$ \bar{U}_N = \sum_{j=1}^{N} U_{N,j}(x, X), $$

where $U_{N,j}$ is exponentially close to $s_j \tanh(x - X_j)/\epsilon$. For $\kappa = 0$ and $\nu = 0$ we have a set of solutions $u$ of (1) that have the form $u = \bar{U}_N(x, X(t)) + \tilde{u}$, where $\tilde{u}$ is a small correction and the kink coordinates $X_j(t)$ evolve in time exponentially slowly: $|dX_j/dt| = O(\exp(-c_0\epsilon^{-1}dist(X)))$ [6]. These solutions are correct while $dist(X) \gg \epsilon$. The time evolution of $X$ is a result of exponentially weak direct kink interaction. For $\kappa \neq 0$ we obtain the kink chain traveling with a constant speed $\kappa: U_N = \bar{U}_N(x - \kappa t, X)$.

In the coming subsection, we consider the case $\nu \neq 0$, where there occurs a non-direct and non-local kink interaction via coupling with $\nu$-reagent.

## 6.2. Equations for kink coordinates

We choose solutions of Eq. (3) having the form $z(x, y, t) = \gamma z_0(x - \kappa t, y)$, where $\gamma > 0$ is a small parameter and $z_0$ is a smooth $2\pi$-periodic function. We substitute $z$ into Eqs. (1), (2) and make variable change $\tilde{x} = x - \kappa t$ (further we omit tilde in notation). Then we obtain the following system

$$ u_t = \frac{\epsilon^2}{2} \Delta u + u - u^3 + \gamma v, \quad (13) $$

$$ v_t = \Delta v + z_0 u_x + \kappa v_x. \quad (14) $$

We are capable to construct asymptotic solutions of the system (13) and (14) under the following assumptions to the small parameters:

**Assumption AP.** Let $\epsilon > 0$ be small enough and

$$ \gamma < \epsilon^4, \quad c_0 \exp(-c_1 \epsilon^{-1/2}) \ll \kappa \ll \gamma, \quad (15) $$

where $c_0, c_1$ are uniform in $\epsilon > 0$.

The main idea in choice of $\gamma$ 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 $\kappa > 0$ should be small as well in order to obtain quasistationary solutions of Eq. (14). The condition $c_0 \exp(-c_1 \epsilon^{-1/2}) \ll \gamma$ is necessary to ensure domination of non-local kink interaction via $\nu$-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 $\gamma > 0$ (see, for example, [6,26]). For sufficiently small $\gamma > 0$ one has

$$ u(x, t) = U_N(x, X(t)) + \tilde{u}(x, t, \gamma), \quad (16) $$

where $U_N(x, X)$ is the kink chain (described above) and $\tilde{u}$ is a correction (the deformation of the kink chain form). The time evolution of $X$ is governed by the equation

$$ \frac{dX_i}{dt} = \gamma (G_i + \tilde{G}_i) + \kappa, \quad (17) $$

where

$$ \tilde{G}_i = (-1)^i \frac{3\tilde{v}_i}{2}, \quad \tilde{v}_i = \int_0^1 v(X_i, y, t) dy \quad (18) $$
and small corrections $\tilde{G}_i$ are uniformly bounded

$$|\tilde{G}_i(X, \gamma, \epsilon)| < c(\gamma^s + \epsilon + \exp(-\epsilon^{-1}\mathrm{dist}(X))),$$

where $s \in (0, 1]$, $c, s$ are uniform in $\gamma$ as $\gamma \to 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

$$\langle f, g \rangle = \int_0^{2\pi} \int_0^1 fg dxdy, \quad \|f\|^2 = \langle f, f \rangle.$$

Let $\Psi_j = \bar{c}_j \epsilon^{-1/2} \cosh^{-2}(\epsilon^{-1}(x - x_i))$, where the normalizing constants $\bar{c}_j$ provide $\|\Psi_j\|_{L_2(\Omega)} = 1$. Note that $\bar{c}_j$ equal $\sqrt{3}/2$ up to exponentially small corrections. The functions $\Psi_j$ are the Goldstone modes induced by kink translation motions.

To obtain the dynamical equations for $X_i$, we impose the condition

$$\langle \tilde{u}(\cdot, \cdot, t), \Psi_j(\cdot - X_i(t)) \rangle = 0 \quad (19)$$

for each $t$ and $i$. These equations define $X$ uniquely for small $\|\tilde{u}\|$, $\gamma > 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 $\tilde{u}$ we then obtain

$$\tilde{u}_t = L_{\epsilon}(X)\tilde{u} + H(\tilde{u}, X, \gamma), \quad (20)$$

where $L_{\epsilon}$ is the linear operator defined by

$$L_{\epsilon}\psi = \frac{\epsilon^2}{2}\Delta\psi + (1 - 3U_N(x, X)^2)\psi, \quad (21)$$

and

$$H(\tilde{u}, X, \gamma) = -\tilde{u}^3 - 3U_N\tilde{u}^2 + \gamma v.$$

The spectrum of the operator $L_{\epsilon}$ is well-studied [6]. Note that $L_{\epsilon}$ 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 $\mathrm{dist}(X) > \delta_1 > 0$ one has

$$\|L_{\epsilon}\Psi_j\| < c_0 \exp(-c_1 \epsilon^{-1}) \quad (22)$$

and if $\langle \psi, \Psi_j \rangle = 0$ for all $j$ then

$$\|L_{\epsilon}\psi\| \le -C_0 \epsilon^2 \|\psi\|, \quad (23)$$

where all constants are uniform in $\epsilon$ as $\epsilon \to 0$.

For a proof see [23].

So, the spectrum of $L_{\epsilon}$ consists of $N$ exponentially small eigenvalues and all the remaining spectrum of $L_{\epsilon}$ lies in the interval $(-\infty, -\delta_0 \epsilon^2)$, where $\delta_0 > 0$ does not depend on $\epsilon > 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_{\epsilon}$ and allows us to solve Eq. (20). It can be done by the standard perturbative methods because for small $\gamma, \kappa$ and $\|\tilde{u}\|$ our model is weakly nonlinear while the linear part is stable due to condition (19).

### 6.3. 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 + \tilde{u}$, where $\tilde{u}$ 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

$$\Delta V_N = -\kappa V_{Nx} - z_0 U_N.$$

We obtain

$$V_N = \sum_{j=1}^{N} W_j(x, y, X), \quad (24)$$

$$\Delta W_j - \kappa W_{jx} = (-1)^j z_0(x, y) U_{N,jx}(x - X_j). \quad (25)$$
To obtain an asymptotic solution it is sufficient to note that for small $\epsilon > 0$ the function $dU_{N,j}/dx = \epsilon^{-1/2}\psi_j + O(\exp(-c\epsilon^{-1}))$ is a good approximation of $\delta$-function (up to a constant uniform in $\epsilon > 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 $\Gamma_m(x)$ the Green function of the one-dimensional boundary value problem satisfying the equation

$$ \frac{d^2 \Gamma_m}{dx^2} - m^2 \Gamma = \delta(x) $$

and the $2\pi$-periodical boundary conditions in x.

Then we resolve (25) by the Fourier method that gives

$$ W_j = \sum_{m=1}^{\infty} \hat{z}_m(X_j) \sin(\pi m y) \Gamma_m(x - X_j), \quad (26) $$

where we remove terms of the order $O(\epsilon)$ and $\hat{z}_m(X_j)$ are the Fourier coefficients of $z_0(X_j, y)$:

$$ \hat{z}_m(X_j) = \frac{b_m}{2} \int_0^1 z_0(X_j, y) \sin(\pi m y) dy, $$

where $m$ are positive integers, $b_m = (1 - (-1)^m)/(\pi m)$ (see Fig. 1).

## 6.4. Hopfield system

Using relations (24), Eqs. (17), (18), and removing small terms, we obtain evolution equations for kink coordinates:

$$ \frac{dX_i}{dt} = \gamma G_i(X), \quad (27) $$

where

$$ G_i(X) = -\frac{2}{3} \sum_{j=1}^{N} \int_{0}^{1} W_j(X_i, y, X_j) dy. \quad (28) $$

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 $\hat{z}_m(X_j)$. First we substitute formula (26) into (28). Then we have

$$ G_i(X) = - \sum_{j=1}^{N} \sum_{m=0}^{\infty} R_{ijm}(X) \quad (29) $$

where

$$ R_{ijm}(X) = \frac{\hat{z}_{2m+1}(X_j) \Gamma_m(X_j - X_i)}{3\pi(2m+1)^2}. $$

The main idea to simplify the formula (29) for $G_i$ is as follows. Suppose that the kinks oscillate at certain fixed points $\bar{X}_j$, i.e.,

$$ X_i(t) = \bar{X}_i + \rho \bar{X}_i(t), \quad (30) $$

where $\bar{X}_i$ are new unknowns and $\rho$ is a small parameter, which is independent of $\epsilon, \gamma, \kappa$ and defines the oscillation magnitude. Suppose temporarily that $|\bar{X}_j|$ stay bounded as $\rho \to 0$ (this assumption will be justified later). We can achieve such behavior of X-solutions of (28) under a special choice of $\hat{z}_{2m+1}(X_j)$. Positions of points $\bar{X}_j$ may be arbitrary but the condition (12) must be satisfied. Namely, for $x$ from a neighborhood of $\bar{X}_j$ we set

$$ \hat{z}_{2m+1}(x) = \xi_{jm} + \eta_{jm}, \quad (31) $$

$$ \xi_{jm}(x) = M_{jm}\sigma(\rho^{-1}(x - \bar{X}_j) - h_j), \quad (32) $$

$$ \eta_{jm}(x) = S_{jm}\rho^{-1}(x - \bar{X}_j), \quad (33) $$

where $\sigma$ is a smooth sigmoidal function, for example, the Fermi function $\sigma(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, $\sigma = \sin$.

We obtain then

$$ R_{ijm}(X) = \bar{R}_{ijm}(X) + \tilde{R}_{ijm}(X), \quad (34) $$
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

where $\tilde{R}_{ijm}(X) = O(\rho)$, are small corrections and

$$R_{ijm}(X) = \Gamma_{ijm}(M_{jm}\sigma(\tilde{X}_j - h_j) + S_{jm}\tilde{X}_j). \quad (35)$$

where $\Gamma_{ijm} = \Gamma_m(\tilde{X}_j - \tilde{X}_i)$

Further, we use the following lemma.

**Lemma II.** For each $N \times N$ matrix $K$ with entries $K_{ij}$, $i = 1, \dots, \bar{N}$ there exist a number $M \ge N$ and coefficients $b_{jm}$ such that

$$\sum_{m=0}^{M} b_{jm} \Gamma_m(\bar{X}_j - \bar{X}_i) = K_{ij} \quad \forall i, j. \quad (36)$$

**Proof.** For the unknown $b_{jm}$ we have a system of linear algebraic equations. For large $m$ we have asymptotics

$$\Gamma_m(x) = (2m)^{-1}(\exp(-m|x|)(1 + o(1))), \quad (37)$$

for $\Gamma_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.
Using this lemma, we can choose $S_{jm}$ and $M_{jm}$ such that $G_i$ take the form

$$G_i(\tilde{X}) = \sum_{j=1}^{N} K_{ij} \sigma(\tilde{X}_j - h_j) - \lambda \tilde{X}_i, \quad (38)$$

where $\lambda > 0$ is a coefficient.

### 6.5. Control of dynamics for Hopfield system

Using (38) we obtain the system (9) for variables $\tilde{X}_i$.

It is easy to show that system (9) has a compact attractor. In fact, $0 \le \sigma(z) \le 1$ thus that system implies the inequalities

$$\frac{dX_i}{dt} \le N|K| - \lambda \tilde{X}_i,$$

where $|K| = \max_{i,j} |K_{ij}|$. These differential inequalities lead to the estimate

$$|\tilde{X}_i(t)| \le (Y_i(0) - N|K|\lambda^{-1}) \exp(-\lambda t) + N|K|\lambda^{-1}.$$

The last estimate shows that system (9) has an absorbing set $A = \{\tilde{X} : |\tilde{X}_i| < N|K|\lambda^{-1}\}$, thus it is dissipative and has a compact attractor. This result justifies our hypothesis on smallness of kink oscillations $\tilde{X}_i$ at points $X_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$, $\lambda$ and $h$.

In coming section, following [13], we show how this result can be applied for cell differentiation.

## 7. 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–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 (rbrrrbbbbr). Suppose that the cell type is determined by a morphogen operator $M(u(\cdot))$ transforming patterns $u(x)$ into patterns $m(x_1), m(x_2), \dots, m(x_n)$, where $x_j$ are cell centers and $m(x_j) \in [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, \dots, s_n)$ consisting of symbols $s_j \in \{r, b\}$ (extension to large alphabets is straightforward). The set of strings $S_D = \{s^{(1)}, s^{(2)}, \dots, s^{(m)}\}$, the time moments $t_j = j\Delta T$ and points $\tilde{x}_1 < \tilde{x}_2 < \dots < \tilde{x}_n$, where $j = 1, \dots, n$ can be called a developmental program if the morphogen values take prescribed values at given time moments in given points:

$$M(u(\cdot, t_k))(x_j) = m(x_j, t_k) = s_j^{(k)}, \quad (39)$$

where $j = 1, \dots, 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.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** This picture illustrates a concept of developmental program. Different types of the cells are shown by different colors. The top row of the cells emerges at the last time moment $t_n$, the previous row appears at $t = t_{n-1}$ and the bottom row arises at the initial time $t = t_0$.</figcaption></figure>

## 8. 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.

### CRediT authorship contribution statement

**Ivan Sudakow:** Conceptualization, Methodology, Writing, Supervision, Visualization, Funding acquisition. **Sergey A. Vakulenko:** Conceptualization, Methodology, Writing, Investigation, Software, Visualization, Funding acquisition. **Dima Grigoriev:** Methodology, Investigation, Formal analysis, Validation.

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

### Acknowledgments

I.S. gratefully acknowledges support from the Division of Physics at the U.S. National Science Foundation (NSF) through Grants PHY-2102906 and in part PHY-1748958. S.V. is supported by the grant of National Institutes of Health (NIH), USA 2R01 OD010936.
# Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cnsns.2022.106844.

## References

[1] Hoffmann LA, Carenza LN, Eckert J, Giomi L. Defect-mediated morphogenesis. 2021, ArXiv:2105.15200v1 [Cond-Mat.Soft].

[2] Keber FC, Loiseau E, Sanchez T, DeCamp JC, Giomi L, Bowick MJ, Marchetti MC, Dogic Z, Bausch AR. Topology and dynamics of active nematic vesicles. Science 2014;345:1135-9.

[3] Livshits A, Shani-Zerbib L, Maroudas-Sacks Y, Braun E, Keren K. Structural inheritance of the actin cytoskeletal organization determines the body axis in regenerating hydra. Cell Reports 2017;18:1410-21.

[4] Luuk Metselaar L, Yeomans JM, Doostmohammadi A. Topology and morphology of self-deforming active shells. Phys Rev Lett 2019;123:208001.

[5] Cheng X, Ferrell JE. Apoptosis propagates through the cytoplasm as trigger waves. 2018, arXiv:1806.04138v1 [q-bio.TO].

[6] Carr J, Pego R. Metastable patterns. Comm on Pure and Appl Math 1989;42:523-76.

[7] Ochta T, Jasnow D. Interface dynamics for layered structures. Phys Rev E 1997;56:5648-58.

[8] Vakulenko S. Dissipative systems generating any structurally stable chaos. Adv Differential Equations 2000;5(7-9):1139-78.

[9] Ruelle D. Elements of differentiable dynamics and bifurcation theory. Elsevier; 2014.

[10] Katok A, Hasselblatt B. Introduction to the modern theory of dynamical systems. Cambridge University Press; 1997.

[11] Moore C. Unpredictability and undecidability in dynamical systems. Phys Rev Lett 1990;64:2354-7.

[12] Moore C. Generalized shifts: Unpredictability and undecidability in dynamical systems. Nonlinearity 1991;4:199-230.

[13] Vakulenko S, Grigoriev D. New way for cell differentiation: reaction, diffusion and chaotic waves. Biosystems 2021. [in press].

[14] FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1961;1:445-66.

[15] Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon. Proc IRE 1962;50:2061-70.

[16] Tao T. On the universality of the incompressible Euler equations on compact manifolds. Discrete Contin Dyn Syst 2018;38(3). http://dx.doi.org/10.3934/dcds.2018064.

[17] Barron AR. Universal approximation bounds for superpositions of a sigmoidal function. Inform Theory, IEEE Transactions on 1993;39(3):930-45.

[18] Henry D. Geometric theory of semilinear parabolic equations. vol. 840, Berlin: Springer-Verlag; 1981.

[19] Constantin P, Foias C, Nicolaenko V, Temam R. Integral manifolds and inertial manifolds for dissipative partial differential equations. In: Springer-verlag, applies mathematical sciences series. 70, Springer; 1988.

[20] Poláčik P. Complicated dynamics in scalar semilinear parabolic equations in higher space dimension. J Differential Equations 1991;89(2):244-71.

[21] Poláčik P. High-dimensional $\omega$-limit sets and chaos in scalar parabolic equations. J Differential Equations 1995;119(1):24-53.

[22] Vakulenko S, Morozov I, Radulescu O. Maximal switchability of centralized networks. Nonlinearity 2016;29(8):2327, URL http://stacks.iop.org/0951-7715/29/i=8/a=2327.

[23] Vakulenko S. Complex attractors and patterns in reaction-diffusion systems. J Dyn Diff Equ 2018;30:175-207.

[24] Kozlov V, Vakulenko S. On chaos in Lotka-Volterra systems: an analytical approach. Nonlinearity 2013;26:2299.

[25] Vakulenko SA. Strange attractors for Oberbeck-Boussinesq model. J Dynam Differential Equations 2021;33(1):303-43.

[26] Vakulenko S, Volpert V. Generalized travelling waves for perturbed monotone reaction-diffusion systems. Nonlinear Analysis. TMA 2001;46:757-76.

[27] https://youtu.be/4frNmeqkhgE, https://youtu.be/4frNmeqkhgE. Video 2021.

[28] Maroto M, Bone RA, J. Kim Dale J. Somitogenesis. Development 2012;139:2453-6.

[29] Resende T, Ferreira M, Teillet MA, Tavares AT, Andrade RP, Palmeirim I. Sonic hedgehog in temporal control of somite formation. Proc Natl Acad Sci USA 2010;107:12907-12.

[30] Pourquié O. The chick embryo: a leading model for model in somitogenesis studies. Mech Dev 2004;121:1069-79.

[31] Cooke J, Zeeman E. A clock and wavefront model for control of the number of repeated structures during animal morphogenesis. J Theoret Biol 1976;58:455-76.

[32] Baker R, Schnell S, Maini P. A clock and wavefront mechanism for somite formation. Dev Biol 2006;293:116-26.