Optimizing mutual synchronization of rhythmic spatiotemporal patterns in reaction-diffusion systems

TL;DR

This paper develops a theoretical framework using phase reduction to optimize the stability of synchronized rhythmic patterns in coupled reaction-diffusion systems, with practical illustrations on FitzHugh-Nagumo models.

Contribution

It derives the optimal linear filter and nonlinear interaction function to maximize synchronization stability in coupled reaction-diffusion systems.

Findings

Derived the optimal linear filter for synchronization stability.

Identified the nonlinear optimal interaction function.

Validated theory with FitzHugh-Nagumo system examples.

Abstract

Optimization of the stability of synchronized states between a pair of symmetrically coupled reaction-diffusion systems exhibiting rhythmic spatiotemporal patterns is studied in the framework of the phase reduction theory. The optimal linear filter that maximizes the linear stability of the in-phase synchronized state is derived for the case where the two systems are linearly coupled. The nonlinear optimal interaction function that theoretically gives the largest linear stability of the in-phase synchronized state is also derived. The theory is illustrated by using typical rhythmic patterns in FitzHugh-Nagumo systems as examples.

Figures

Equations

$$\frac{\partial }{\partial t}{\mathbf{X}}_1(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_2(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_1(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_2(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_1(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_2(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_1(\mathbf{r},t)$$

$$\frac{\partial }{\partial t}{\mathbf{X}}_2(\mathbf{r},t)$$

$$\frac{\partial }{\partial \theta }\mathbf{Z}(\mathbf{r},\theta )=-D\mathbf{F}({\mathbf{X}}_0(\mathbf{r},\theta ),\mathbf{r}{)}^{†}\mathbf{Z}(\mathbf{r},\theta )-{\hat{D}}^{†}{\nabla }^2\mathbf{Z}(\mathbf{r},\theta )$$

$${\int }_{V}d\mathbf{r}\mathbf{Z}(\mathbf{r},\theta )\cdot \mathbf{U}(\mathbf{r},\theta )=1$$

$${\dot{\theta }}_1(t)$$

$${\dot{\theta }}_2(t)$$

$$\Gamma (\varphi )=\frac{1}{2\pi }{\int }_0^{2\pi }d\psi {\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\mathbf{Z}\left(\mathbf{r},\psi +\varphi \right)\cdot \hat{A}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\mathbf{X}}_0\left({\mathbf{r}}^{\prime },\psi \right),$$

$$\Gamma (\varphi )=\frac{1}{2\pi }{\int }_0^{2\pi }d\psi {\int }_{V}d\mathbf{r}\mathbf{Z}\left(\mathbf{r},\psi +\varphi \right)\cdot {\mathbf{X}}_0\left(\mathbf{r},\psi \right).$$

$$\dot{\varphi }(t)=ϵ{\Gamma }_a(\varphi ),$$

$${\Gamma }_a(\varphi )=\Gamma (\varphi )-\Gamma (-\varphi )$$

$${\Gamma }_a^{\prime }(0)={\frac{d}{d\varphi }{\Gamma }_a(\varphi )|}_{\varphi =0}.$$

$${\Gamma }_a^{\prime }(0)=2{\Gamma }^{\prime }(0)$$

$$=-2\cdot \frac{1}{2\pi }{\int }_0^{2\pi }d\psi {\int }_{V}d\mathbf{r}\mathbf{Z}\left(\mathbf{r},\psi \right)\cdot \mathbf{U}\left(\mathbf{r},\psi \right)=-2,$$

$${\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }{\Vert \hat{A}(\mathbf{r},{\mathbf{r}}^{\prime })\Vert }_F^2=P,$$

$${\Vert \hat{A}(\mathbf{r},{\mathbf{r}}^{\prime })\Vert }_F^2=\sum _{i=1}^m\sum _{j=1}^m{A}_{ij}(\mathbf{r},{\mathbf{r}}^{\prime }{)}^2.$$

$$S\{\hat{A},\lambda \}=-{\Gamma }_a^{\prime }(0)-\lambda \left({\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{A}_{ij}{\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}^2-P\right),$$

$${\Gamma }^{\prime }(\varphi )=\frac{1}{2\pi }{\int }_0^{2\pi }d\psi {\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\frac{\partial \mathbf{Z}\left(\mathbf{r},\psi +\varphi \right)}{\partial \psi }\cdot \hat{A}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\mathbf{X}}_0\left({\mathbf{r}}^{\prime },\psi \right).$$

$${\Gamma }_a^{\prime }(0)=2{\Gamma }^{\prime }(0)$$

$$=-2\cdot \frac{1}{2\pi }{\int }_0^{2\pi }d\psi {\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\mathbf{Z}\left(\mathbf{r},\psi \right)\cdot \hat{A}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)\mathbf{U}\left({\mathbf{r}}^{\prime },\psi \right).$$

$${\Gamma }_a^{\prime }(0)$$

$$=-2{\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{A}_{ij}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)W_{ij}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right).$$

$$W_{ij}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)=\frac{1}{2\pi }{\int }_0^{2\pi }d\psi Z_i\left(\mathbf{r},\psi \right)U_j\left({\mathbf{r}}^{\prime },\psi \right),$$

$$S\{\hat{A},\lambda \}$$

$$-\lambda \left({\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{A}_{ij}{\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}^2-P\right),$$

$$W_{ij}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)=\lambda A_{ij}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right),$$

$$\hat{A}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{\lambda }\hat{W}(\mathbf{r},{\mathbf{r}}^{\prime }).$$

$${\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{A}_{ij}{\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}^2=P.$$

$$\lambda =\sqrt{\frac{1}{P}{\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{W}_{ij}{\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}^2},$$

$${\Gamma }_a^{\prime }(0)=-\frac{2}{\lambda }{\int }_{V}d\mathbf{r}{\int }_{V}d{\mathbf{r}}^{\prime }\sum _{i=1}^m\sum _{j=1}^m{W}_{ij}{\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)}^2.$$

Full text

I Introduction

Synchronization of rhythmic systems is widely observed in the real world and has been extensively studied in various areas of science and engineering ref:winfree80 ; Kuramoto ; Glass ; ref:pikovsky01 ; ref:hoppensteadt97 ; ref:ermentrout10 ; bullo . In biological systems, synchronization often plays significant functional roles, such as generation of heartbeats, rhythmic gaits, and circadian rhythms. In some engineering systems, synchronization is a precondition for their functionality, where power grids provide a well-known example.

Recently, synchronization between non-conventional self-oscillatory systems has recently attracted considerable attention, for example, dynamical systems with time-delayed feedback Kotani ; Novicenko , spatially extended reaction-diffusion systems Fukushima ; Hildebrand ; Nakao , and fluid systems read2010 ; read2012 ; Kawamura2 . In analyzing synchronization properties of rhythmic systems described as weakly perturbed limit-cycle oscillators, the phase reduction theory has been used as a standard method for clarifying their mechanisms  Kuramoto ; ref:winfree80 ; Glass ; ref:pikovsky01 ; ref:hoppensteadt97 ; ref:ermentrout10 ; ref:brown04 ; ref:nakao16 ; ref:ashwin16 . Recent developments in the phase reduction theory Kotani ; Novicenko ; Nakao ; Kawamura2 have shown that, even if the state space of the rhythmic system is infinite-dimensional, it can still exhibit synchronization in a similar way to low-dimensional oscillators as long as it exhibits stable limit-cycle oscillation.

In this study, we consider synchronization of reaction-diffusion systems exhibiting rhythmic spatiotemporal patterns. Reaction-diffusion systems have played important roles in modeling a variety of spatiotemporal patterns that arise in chemical and biological systems ref:mikhailov06 ; ref:mikhailov13 ; ref:mikhailov14 ; ref:cross93 ; ref:cross09 ; ref:engel14 ; ref:engel16 . Among them, rhythmic spatiotemporal patterns, such as oscillating spots, target waves, and rotating spirals, can be regarded as stable limit-cycle oscillations of the reaction-diffusion systems. Synchronization between rhythmic spatiotemporal patterns has been experimentally realized using coupled electrochemical systems exhibiting reaction waves of ${\rm H}_{2}{\rm O}_{2}$ reduction on ${\rm Pt}$ ring electrodes, where two waves are coupled via the common chemical solution Fukushima , and coupled photosensitive Belousov-Zhabotinsky systems exhibiting spiral patterns, where the two patterns are coupled via video cameras and projectors Hildebrand (see also Sakurai ; ref:mikhailov06 ; Epstein ; ref:showalter15 ; ref:weiss15 ; ref:weiss17 ). Synchronization of rhythmic fluid flows has also been studied and its possible importance in global climate has been argued read2010 ; read2012 .

In our recent work Nakao , we generalized the conventional phase reduction theory for finite-dimensional limit-cycle oscillators to limit-cycle oscillations of reaction-diffusion systems with infinite-dimensional state space. Using the theory, we derived the phase sensitivity function, which characterizes linear phase response of the rhythmic pattern to weak perturbations, and analyzed mutual synchronization between a pair of reaction-diffusion systems coupled by linear diffusive interaction. We also developed similar phase reduction theories for the collective oscillations in globally coupled noisy oscillators Kawamura1 and for oscillatory thermal convection in a Hele-Shaw cell Kawamura2 . Moreover, we analyzed synchronization between non-interacting convection cells exhibiting oscillatory thermal convection caused by common noise and derived the optimal input pattern for stable noise-induced synchronization Kawamura3 .

In this study, we consider the case that a pair of reaction-diffusion systems, both of which are exhibiting rhythmic spatiotemporal patterns, are mutually coupled via weak symmetric interaction. In the case of the simplest linear diffusive interaction where every point of the system is coupled to a corresponding point of the other system, we have shown that the two systems can undergo mutual synchronization Nakao . However, the phase sensitivity function of the rhythmic pattern is often strongly localized in space, and introducing interaction at every point in the system as in Ref. Nakao would not be generally efficient. In this study, we aim to clarify the theoretical limit of efficiency in synchronizing rhythmic patterns via mutual coupling by seeking for optimal interaction schemes that realize stable in-phase synchronization between two reaction-diffusion systems in the framework of the phase reduction theory.

Regarding optimization of synchronization, optimal input signals that efficiently entrains a limit-cycle oscillator described by ordinary differential equations have been obtained for various situations Zlotnik ; ref:moehlis06 ; ref:harada10 ; ref:dasanayake11 ; ref:zlotnik12 ; ref:zlotnik12 ; ref:pikovsky15-prl ; ref:tanaka14a ; ref:tanaka14b ; ref:tanaka15 ; ref:zlotnik2016 . Also, in our preceding article, we derived optimal cross-coupling matrices that maximize linear stability of the synchronized states in a pair of diffusively coupled limit-cycle oscillators described by ordinary differential equations Shirasaka . In this paper, we further generalize the analysis to a pair of coupled reaction-diffusion systems exhibiting rhythmic patterns and try to derive the optimal interaction function. We first restrict ourselves to a practical situation where the interaction between the two systems is linear and derive the optimal filtering function for stable synchronization. We then derive the optimal nonlinear interaction function between the systems to clarify the theoretical limit to the improvement of stability. The results are illustrated by using rhythmic spatiotemporal patterns in FitzHugh-Nagumo reaction-diffusion systems, that is, a traveling pulse and an oscillating spot in one dimension, and a rotating spiral in two dimensions.

II Theory

II.1 A pair of mutually coupled reaction-diffusion systems

We consider a pair of mutually coupled $m$-component reaction-diffusion systems in a $d$-dimensional space exhibiting stable limit-cycle oscillations (see Figures for typical rhythmic spatiotemporal patterns of the FitzHugh-Nagumo reaction-diffusion systems), described by

[TABLE]

in some spatial domain $V\subset{\bm{R}}^{d}$. Here, ${\bm{r}}\in{\bm{R}}^{d}$ is the spatial location, $t\in{\bm{R}}$ is the time, ${\bm{X}}_{1,2}:{\bm{R}}^{d}\times{\bm{R}}\to{\bm{R}}^{m}$ are the spatial patterns of the systems, i.e., the system states, ${\bm{F}}:{\bm{R}}^{m}\times{\bm{R}}^{d}\to{\bm{R}}^{m}$ represents the dynamics of the system, $\hat{D}\in{\bm{R}}^{m\times m}$ is a matrix of diffusion constants, and $\nabla^{2}$ is the Laplacian operator. For simplicity, we consider two identical systems whose dynamics are described by the same function ${\bm{F}}$ and diffusion matrix $\hat{D}$, and also assume that they are symmetrically coupled. The dynamics ${\bm{F}}$ can depend on the location ${\bm{r}}$, for example, the excitability of the system can be different from place to place. The last term in the right-hand side of each equation represents mutual linear interaction between the two systems, where each system is coupled to the other system via a $m\times m$ matrix of spatial linear filters ${\hat{A}}({\bm{r}},{\bm{r}}^{\prime}):{\bm{R}}^{d}\times{\bm{R}}^{d}\to{\bm{R}}^{m\times m}$. The parameter $\epsilon\geq 0$ is the interaction intensity, which is assumed to be small.

As a benchmark, we also consider the following systems with simple, direct mutual interaction:

[TABLE]

where every point in the system is directly coupled to the corresponding point of the other system without filtering. The definitions of the variables and parameters are the same as in Eq. (2). Both interaction schemes, Eq. (2) and Eq. (4), can exhibit in-phase synchronization between the systems, and we compare the stability of the synchronized states between them.

As another, typical interaction scheme, we may also consider diffusive interaction between the two reaction-diffusion systems given by

[TABLE]

It is clear that this interaction scheme also allows in-phase synchronization. For sufficiently small $\epsilon$, which is assumed throughout this study, we can show that the stability of the in-phase synchronized state with this diffusive interaction scheme is approximately equal (up to $O(\epsilon)$) to that for the interaction scheme given by Eq. (2). Similarly, in our previous paper Nakao , we analyzed the following simple case with direct diffusive interaction, where every point in the system is diffusively coupled to the corresponding point in the other system as

[TABLE]

and showed that the two systems undergo mutual synchronization by using the phase reduction theory. Linear stability of the in-phase synchronized state of this interaction scheme is also approximately the same as that for Eq. (4) when $\epsilon$ is sufficiently small.

In this study, we focus on the interaction schemes given by Eq. (2) and Eq. (4) and analyze their synchronization properties. We consider a general nonlocal interaction given by Eq. (2) and try to optimize the linear filter $\hat{A}({\bm{r}},{\bm{r}}^{\prime})$ so that the two systems exhibit more stable in-phase synchronization than the case with the simple interaction given by Eq. (4). In the following, we refer to the interaction scheme in Eq. (2) as nonlocal, while that in Eq. (4) as direct.

II.2 Phase reduction

In Ref. Nakao , we generalized the phase reduction theory for finite-dimensional limit-cycle oscillators Kuramoto to reaction-diffusion systems exhibiting stable rhythmic patterns. Using the theory, we can systematically approximate the dynamics of two weakly coupled reaction-diffusion systems by simple two-dimensional coupled phase equations and analyze their synchronization properties.

Suppose that a single reaction-diffusion system (without interaction) exhibits a stable rhythmic pattern, that is, a stable limit-cycle solution ${\bm{X}}_{0}({\bm{r}},t+T)={\bm{X}}_{0}({\bm{r}},t)$ of period $T$ and frequency $\omega=2\pi/T$. We can introduce a phase variable $\theta\in[0,2\pi]$ of the rhythmic pattern around the limit-cycle solution in the state space of the system that always increases with a constant frequency $\omega$ in the absence of perturbation (i.e., the asymptotic phase ref:winfree80 ; Kuramoto ; ref:hoppensteadt97 ; ref:ermentrout10 ; ref:brown04 ; ref:nakao16 ; ref:ashwin16 ), and represent the limit-cycle solution as a function of phase $\theta$, rather than time $t$, as ${\bm{X}}_{0}({\bm{r}},\theta)$ ($0\leq\theta\leq 2\pi$).

In the phase reduction theory, the phase sensitivity function ${\bm{Z}}({\bm{r}},\theta):{\bm{R}}^{d}\times[0,2\pi]\to{\bm{R}}^{m}$ plays an important role, which characterizes linear phase response of the system to a weak perturbation that is applied when the system state is at ${\bm{X}}_{0}({\bm{r}},\theta)$. It is given by a $2\pi$-periodic solution (an eigenfunction associated with the zero eigenvalue) to the adjoint equation

[TABLE]

with appropriate boundary conditions, where $D{\bm{F}}$ is the Jacobi matrix of ${\bm{F}}({\bm{X}},{\bm{r}})$ at ${\bm{X}}={\bm{X}}_{0}({\bm{r}},\theta)$ and ${\dagger}$ denotes matrix transpose, and satisfies a normalization condition

[TABLE]

for $0\leq\theta\leq 2\pi$. Here, we defined the tangent field ${\bm{U}}({\bm{r}},\theta)$ of ${\bm{X}}_{0}({\bm{r}},\theta)$ along the limit-cycle solution as ${\bm{U}}({\bm{r}},\theta)={\partial{\bm{X}}_{0}({\bm{r}},\theta)}/{\partial\theta}$.

We consider weakly coupled reaction-diffusion systems given by Eq. (2) or Eq. (4), and assume that the rhythmic patterns are only slightly perturbed and persist even when weak mutual interaction between the two systems is introduced. We can then approximately describe the system states using only scalar phase variables $\theta_{1,2}\in[0,2\pi]$ as ${\bm{X}}_{1,2}({\bm{r}},t)={\bm{X}}_{0}({\bm{r}},\theta_{1,2}(t))$, and derive approximate phase equations for $\theta_{1,2}(t)$ from Eq. (2) or Eq. (4) as

[TABLE]

where the overdot ($\dot{}$) represents $d/dt$ and the $2\pi$-periodic function $\Gamma(\phi):[0,2\pi]\to{\bm{R}}$ is called the phase coupling function. In the case of the nonlocal interaction, Eq. (2), the phase coupling function is given by

[TABLE]

and in the case of the direct interaction, Eq. (4), the phase coupling function is simply given by

[TABLE]

Synchronization between the two systems can be analyzed in the same way as for finite-dimensional coupled oscillators Kuramoto . From Eq. (12), we can derive the dynamics of the phase difference $\phi=\theta_{1}-\theta_{2}$, which we restrict in the range $[-\pi,\pi]$, as

[TABLE]

where

[TABLE]

is the antisymmetric part of the phase coupling function $\Gamma(\phi)$. Fixed points of this equation satisfying $\Gamma_{a}(\phi^{*})=0$ correspond to the phase differences where the two systems exhibit synchronization. The stability of a fixed point $\phi=\phi_{*}$ is characterized by the slope of $\Gamma_{a}(\phi_{*})$.

Because the function $\Gamma_{a}(\phi)$ vanishes at $\phi=0$ and $\phi=\pm\pi$ by definition, Eq. (15) has fixed points at $\phi^{*}=0$ and $\phi^{*}=\pm\pi$. We focus on the in-phase synchronized state, $\phi^{*}=0$, whose stability is characterized by

[TABLE]

In the case of the nonlocal interaction, Eq. (2), this value depends on ${\hat{A}}({\bm{r}},{\bm{r}}^{\prime})$. For the direct interaction, Eq. (4), this value can be explicitly calculated as

[TABLE]

where we performed partial integration, used the $2\pi$-periodicity of ${\bm{Z}}$ and ${\bm{X}}_{0}$ to eliminate the surface terms, and used the normalization condition Eq. (10) for ${\bm{Z}}$. In this study, we try to make $\Gamma_{a}^{\prime}(0)$ as negative as possible in order to improve the stability of the in-phase synchronization under a constraint for the norm of ${\hat{A}}({\bm{r}},{\bm{r}}^{\prime})$.

II.3 Optimal linear filter for stable synchronization

Under the framework of the phase reduction approximation for sufficiently small $\epsilon$, we seek for the optimal filter function ${\hat{A}}({\bm{r}},{\bm{r}}^{\prime})$ that maximizes the linear stability $-\Gamma^{\prime}_{a}(0)$ of the in-phase synchronized state. As a constraint, we fix the spatial average of the Frobenius norm (or the Hilbert-Schmidt norm) of the linear filter $A({\bm{r}},{\bm{r}}^{\prime})$, i.e.,

[TABLE]

where $P>0$ is a constant and

[TABLE]

To this end, we consider the action

[TABLE]

where $\lambda$ is the Lagrange multiplier, and find $\hat{A}$ and $\lambda$ that give the extremum of $S$.

From Eq. (13), the slope $\Gamma^{\prime}(\phi)=d\Gamma(\phi)/d\phi$ of the phase coupling function is given by

[TABLE]

Thus, the slope $\Gamma_{a}^{\prime}(0)$ of the antisymmetric part $\Gamma_{a}(\phi)=\Gamma(\phi)-\Gamma(-\phi)$ at $\phi=0$ is calculated as

[TABLE]

Denoting the vector components of ${\bm{Z}}$, ${\bm{X}}_{0}$, and ${\bm{U}}=\partial{\bm{X}_{0}}/\partial\psi$ as $Z_{1},\cdots,Z_{m}$, $X_{1},\cdots,X_{m}$, and $U_{1},\cdots,U_{m}$, respectively, and the matrix components of $\hat{A}$ as $\{A_{ij}\}$ ($i,j=1,2,\cdots,m$), $\Gamma_{a}^{\prime}(0)$ can be expressed as

[TABLE]

Here, we have defined a correlation matrix $\hat{W}\left(\bm{r},\bm{r}^{\prime}\right)$ whose components are given by

[TABLE]

which characterizes the spatial correlation between ${\bm{Z}}$ and ${\bm{U}}$ averaged over one period of oscillation.

The action is now given by

[TABLE]

and by taking variations with respect to $A_{ij}({\bm{r}},{\bm{r}}^{\prime})$, we obtain

[TABLE]

so the optimal linear filter is given by

[TABLE]

Differentiating $S$ with respect to $\lambda$ simply gives the constraint,

[TABLE]

This equation gives the Lagrange multiplier $\lambda$ as

[TABLE]

where we have chosen the plus sign so that $S$ takes the extremum (or $\Gamma^{\prime}(0)$ becomes negative) at the optimal $\hat{A}$. The largest negative slope of $\Gamma_{a}^{\prime}(0)$ is given by

[TABLE]

This value gives the largest stability of the in-phase synchronized state for the nonlocal interaction.

Note that spatial linear filtering of the field ${\bm{X}}$ by $\hat{A}$ in Eq. (2) is a fundamental method of image processing and can be easily performed. It is also notable that the expression for the correlation matrix, Eq. (28), consists of only two terms, ${\bm{U}}$ and ${\bm{Z}}$. These are the most fundamental quantities of limit-cycling systems, that is, ${\bm{U}}$ is a tangent field to the limit-cycle orbit and ${\bm{Z}}$ is a phase sensitivity function of a limit-cycle orbit, which are adjoint to each other. These quantities always arise in the phase reduction analysis of coupled oscillators and the expression of Eq. (28) is physically natural.

II.4 Optimal nonlinear interaction

Theoretically, we can consider a more general case that the two system are coupled via nonlinear functional as

[TABLE]

where ${\bm{H}}:C\times{\bm{R}}^{d}\to{\bm{R}}^{m}$ ($C$ represents the set of spatial patterns) is a functional of the spatial pattern. The phase coupling function $\Gamma(\phi):[0,2\pi]\to{\bm{R}}$ in this case can be calculated as

[TABLE]

where ${\bm{H}}({\bm{r}},\psi)={\bm{H}}\{{\bm{X}}_{0}(\cdot,\psi),{\bm{r}}\}$. Note that we have replaced the functional ${\bm{H}}\{{\bm{X}}_{0}(\cdot,\psi),{\bm{r}}\}$ of ${\bm{X}}_{0}({\bm{r}},\psi)$ by the function ${\bm{H}}({\bm{r}},\psi)$ of $\psi$ here, because the function ${\bm{X}}_{0}({\bm{r}},\psi)$ is solely determined by phase $\psi$ in the framework of the phase reduction. We refer to ${\bm{H}}({\bm{r}},\psi)$ as the interaction function.

We try to make $\Gamma_{a}^{\prime}(0)$ as negative as possible under the constraint

[TABLE]

that is, we fix the squared mean of ${\bm{H}}$ over the space and all possible combinations of the phase variables at some constant $Q>0$. This amounts to fixing the average “energy” of the mutual interaction between the two systems. We consider the following action functional of ${\bm{H}}$ and a Lagrange multiplier $\lambda$:

[TABLE]

where $\Gamma_{a}^{\prime}(0)$ in the first term can be represented as

[TABLE]

By taking variations of $S\{{\bm{H}},\lambda\}$ with respect to ${\bm{H}}$, we obtain a Euler-Lagrange equation

[TABLE]

which yields

[TABLE]

The Lagrange multiplier $\lambda$ is given by

[TABLE]

Note that the plus sign, which gives the extremum of $S\{\bm{H},\lambda\}$, has been chosen here. The largest negative slope of the antisymmetric part of the phase coupling function at $\phi=0$ is given by

[TABLE]

which yields the largest possible stability of the in-phase synchronized state for general nonlinear interaction of Eq. (37).

Thus, ${\bm{H}}({\bm{r}},\psi)\propto-\partial{\bm{Z}}({\bm{r}},\psi)/\partial\psi$ is the optimal interaction function in the nonlinear case. That is, the optimal nonlinear interaction between the two systems is realized by (i) measuring the phase $\psi$ of the other system and (ii) driving the system using the negative derivative of the phase sensitivity function with respect to phase $\psi$. This result is consistent with Zlotnik et al.’s Zlotnik for the optimal periodic input signal that maximizes linear stability of the entrainment of ordinary limit-cycle oscillators.

In practice, however, online continuous-time estimation of the instantaneous phase value $\psi$ from an observed rhythmic pattern is generally not straightforward (note that we need to estimate the correct asymptotic phase of a given pattern). For low-dimensional oscillators, we could construct a mapping from the oscillator states to phase values beforehand and use it to estimate the instantaneous phase value of the observed oscillator state in real time. However, construction of such a mapping from high-dimensional data of spatial patterns to phase values is generally a non-trivial, demanding task. Thus, the optimal nonlinear interaction may not be easy to realize experimentally, in contrast to the the nonlocal interaction that can easily be implemented once the optimal linear filter $\hat{A}$ is given.

III Numerical simulations

III.1 FitzHugh-Nagumo model

We now illustrate the theoretical results with numerical examples. As the reaction-diffusion system, we use the FitzHugh-Nagumo (FHN) system described by two field variables, ${\bm{X}}=(X_{u},X_{v})=(u,v)$, which obey

[TABLE]

This system can exhibit various rhythmic spatiotemporal patterns such as traveling pulses, oscillating spots, target waves, and rotating spirals. Here, as typical examples, we consider a traveling pulse (on a ring) and an oscillating spot in one-dimensional systems, and a rotating spiral in two-dimensional systems. The setup of simulations are basically the same as in Ref. Nakao , but some of the parameter values and system sizes are modified.

We compare the results for the nonlocal interaction, Eq. (2) with the optimal linear filter $\hat{A}$, the direct interaction, Eq. (4), and the optimal nonlinear interaction, Eq. (37) with the optimal function ${\bm{H}}$. To make a fair comparison between different interaction schemes, we fix the squared mean of the interaction term over the spatial domain and over the phase to a constant, namely, we calculate the quantity

[TABLE]

for the direct interaction and appropriately normalize the nonlocal interaction,

[TABLE]

and the nonlinear interaction,

[TABLE]

so that

[TABLE]

and

[TABLE]

are satisfied. That is, we fix the average “energy” of the interaction functions between the two systems over one period of oscillation.

For each pattern, the limit-cycle solution, phase sensitivity function, and antisymmetric part of the phase coupling functions for the direct, optimal nonlocal, and optimal nonlinear interactions are shown. Direct numerical simulations of the synchronization process of the reaction-diffusion systems are also shown for the direct and optimal nonlocal interactions (see Appendix on numerical computation of the optimal nonlocal interaction). Optimal nonlinear interaction is not simulated, because it requires the instantaneous phase values of the spatial patterns and therefore is not easy to realize. We use the synchronization error, $E=\sqrt{\int_{V}d{\bm{r}}\,\|{\bm{X}}_{1}({\bm{r}},t)-{\bm{X}}_{2}({\bm{r}},t)\|^{2}}$, and the phase difference, $\phi=\theta_{1}-\theta_{2}$, measured from the simulated patterns to see the convergence of the systems to synchronization. Here, the phase difference is measured stroboscopically by using threshold-crossing times of the patterns roughly at intervals of the period of oscillation.

III.2 Traveling pulse

Figures 1-4 show the results for a traveling pulse in a one-dimensional space $0\leq x\leq L$ with periodic boundary conditions. The parameters are $\alpha=0.1$, $\tau^{-1}=0.002$, $\gamma=2.5$, $D_{u}=1.0$, and $D_{v}=0.1$. The system size is $L=200$ and discretized by using $N=200$ grid points. The interaction intensity between the systems is $\epsilon=0.0001$.

Figure 1 shows the traveling-pulse solution ${\bm{X}}_{0}(x,\theta)$, the tangent function ${\bm{U}}(x,\theta)=\partial{\bm{X}}_{0}(x,\theta)/\partial\theta$, and the phase sensitivity function ${\bm{Z}}(x,\theta)$. In each figure, both $u$ and $v$ components of ${\bm{X}}_{0}(x,\theta=0)=(X_{u},X_{v})$, ${\bm{U}}(x,\theta=0)=(U_{u},U_{v})$, and ${\bm{Z}}(x,\theta=0)=(Z_{u},Z_{v})$ are plotted. Because the pulse keeps a constant shape and simply travels to the right with a constant velocity, all the functions simply translate to the right with a constant velocity without changing their shapes. The period of the oscillation is $T\approx 395.7$ and the frequency is $\omega\approx 0.01588$. From these data, the correlation matrix $\hat{W}(x,x^{\prime})$ is calculated.

Figure 2 shows the nonlocal interaction ${\bm{G}}(x,\theta)$ and the nonlinear interaction ${\bm{H}}(x,\theta)$ at $\theta=0$. These optimal interaction functions also translate to the right with the pulse without changing their shapes. The nonlocal interaction function and nonlinear interaction function differ from each other. However, they have one thing in common; that is, both interaction functions change their signs in front of the pulse. This is actually essential for efficient control of the phase of the pulse.

Figure 3 shows the antisymmetric part $\Gamma_{a}(\phi)$ of the phase coupling function for the direct, nonlocal optimal, and optimal nonlinear interaction. Both nonlocal interaction and nonlinear interaction yield much higher stability of the synchronized state than the direct interaction. The nonlinear interaction gives the highest linear stability, but the nonlocal interaction also yields reasonably high stability. We can also observe that both nonlocal interaction and nonlinear interaction uniquely give the in-phase synchronized state as the globally stable solution, while the direct interaction yields both in-phase and anti-phase synchronized solutions as stable solutions. Because global stability of the solutions is not considered in the present optimization, these results are coincidental.

Figure 4 shows the synchronization dynamics between the two reaction-diffusion systems obtained by direct numerical simulations for the direct and optimal nonlocal interactions. Temporal evolution of the synchronization error and the phase difference (measured stroboscopically at each period of oscillation) are shown. We can clearly see that the optimal nonlocal interaction yields much faster synchronization than the direct coupling. When the phase difference is sufficiently small, exponential growth (or decay) rate of the phase difference coincides with the slope $\Gamma_{a}^{\prime}(0)$.

III.3 Oscillating spot

Figures 5-10 show the results for the oscillating-spot solution in a one-dimensional space $[0,L]$ with no-flux boundary conditions. In this case, the parameter $\alpha(x)$ that controls excitability of the media is spatially modulated as $\alpha(x)=\alpha_{0}+(\alpha_{1}-\alpha_{0})(2x/L-1)^{2}$ with $\alpha_{0}=-1.1$ and $\alpha_{1}=-1.6$, so that the oscillating spot stays at the center of the system. The other parameters are $\tau^{-1}=0.028$, $\gamma=2.0$, $D_{u}=1.0$, and $D_{v}=2.5$. The length of the system is $L=80$ and discretized by using $N=240$ grid points. The intensity of mutual interaction between the systems is $\epsilon=0.00001$.

Figure 5 shows the oscillating-spot solution ${\bm{X}}_{0}(x,\theta)$, the tangent function ${\bm{U}}(x,\theta)=\partial{\bm{X}}_{0}(x,\theta)/\partial\theta$, and the phase sensitivity function ${\bm{Z}}(x,\theta)$, all at $\theta=0$. In each figure, both $u$ and $v$ components of ${\bm{X}}_{0}(x,\theta=0)=(X_{u},X_{v})$, ${\bm{U}}(x,\theta=0)=(U_{u},U_{v})$, and ${\bm{Z}}(x,\theta=0)=(Z_{u},Z_{v})$ are plotted. The period of the oscillation is $T\approx 200.6$ and the frequency is $\omega\approx 0.0313$. Figure 6 shows the evolution of ${\bm{X}}_{0}(x,\theta)$ and ${\bm{Z}}(x,\theta)$ in color code for one period of oscillation, $0\leq\theta\leq 2\pi$. The phase sensitivity function is strongly localized near the interfaces of the spot.

Figure 7 shows all four components of the correlation matrix $\hat{W}(x,x^{\prime})$, and Fig. 8 shows the optimal nonlocal interaction ${\bm{G}}(x,\theta)$ and the optimal nonlinear interaction ${\bm{H}}(x,\theta)$ for one period of oscillation $(0\leq\theta\leq 2\pi)$. The results for the nonlocal and nonlinear cases are different, but the locations at which the interaction functions change their signs are roughly similar in both cases and reflect the locations where the phase sensitivity functions take large values.

Figure 9 shows the antisymmetric part $\Gamma_{a}(\phi)$ of the phase coupling function for the direct, optimal nonlocal, and optimal nonlinear interaction. Both nonlocal interaction and nonlinear interaction yield much higher stability than the direct interaction. In this case, the optimal nonlocal interaction yields somewhat lower stability than the optimal nonlinear case. This discrepancy arises because, in addition to the inevitable smoothness of the interaction function due to filtering in the nonlocal case, the spatial linear filter $\hat{A}$ cannot shift the temporal phase of the interaction function ${\bm{G}}$ from the phase of the oscillation of the spot, while the phase of the optimal nonlinear interaction ${\bm{H}}$ is slightly shifted from that of the spot. In this example, while in-phase and anti-phase synchronized solutions are both stable for the direct interaction, both nonlocal interaction and nonlinear interaction give global stability of the in-phase synchronized solution as in the previous case of traveling pulses. Figure 10 shows the evolution of the synchronization error and the phase difference between the systems for the direct and optimal nonlocal interactions. Optimal nonlocal interaction yields much faster convergence to the synchronized state.

III.4 Rotating spiral

This example is motivated by the experimental study by Hildebrand et al. Hildebrand . Figures 11-14 show the results for the spiral in a two-dimensional square $0\leq x,y\leq L$ with no-flux boundary conditions. The parameter $\alpha(x)$ that controls excitability of the media is spatially modulated as $\alpha(x,y)=\alpha_{0}+(\alpha_{1}-\alpha_{0})\exp(-r^{4}/r_{0}^{4}),r=\sqrt{(x-L/2)^{2}+(y-L/2)^{2}}$ with $\alpha_{0}=0.05$ and $\alpha_{1}=0.5$, so that the spiral stays at the center of the system. The other parameters are $\tau^{-1}=0.005$, $\gamma=2.5$, $D_{u}=1.0$, and $D_{v}=0.0$. The size of the system is $L\times L=80\times 80$ and discretized by using $N^{2}=80^{2}$ grid points. The interaction intensity between the systems is $\epsilon=0.00001$.

Figure 11 shows the spiral solution ${\bm{X}}_{0}(x,y,\theta)$ and the phase sensitivity function ${\bm{Z}}(x,y,\theta)$ at $\theta=0$. The spiral pattern ${\bm{X}}_{0}(x,y,\theta)$ keeps a constant shape and rotates around the center with a constant frequency and, accordingly, ${\bm{U}}(x,y,\theta)$ (not shown) and ${\bm{Z}}(x,y,\theta)$ also rotates around the center. The period of the oscillation is $T\approx 179$ and the frequency is $\omega\approx 0.0351$. It can be seen that the phase sensitivity function ${\bm{Z}}(x,y,\theta)$ is strongly localized near the spiral tip.

Figure 12 shows the phase sensitivity function ${\bm{Z}}(x,y,\theta)$, optimal nonlocal interaction function ${\bm{G}}(x,y,\theta)$, and the optimal nonlinear interaction function ${\bm{H}}(x,y,\theta)$, enlarged near the spiral tip at $\theta=0$. The patterns of ${\bm{G}}$ and ${\bm{H}}$ are different, but they are similar in that both of them are localized near the spiral tip and exhibit localized positive and negative spots.

Figure 13 shows the antisymmetric part $\Gamma_{a}(\phi)$ of the phase coupling function for the direct, optimal nonlocal, and optimal nonlinear interaction. Both nonlocal interaction and nonlinear interaction give drastically higher stability than the direct interaction. The optimal nonlinear interaction yields the highest stability, but the optimal nonlocal interaction also yields reasonably high stability. In this example, again, the nonlocal interaction and the optimal nonlinear interaction give global stability of the in-phase synchronized solution, while the direct interaction provides bistability of in-phase and anti-phase synchronized solutions.

Finally, Fig. 14 shows the evolution of the synchronization error and the phase difference between the systems for the direct and optimal nonlocal interactions. Optimal nonlocal interaction yields much faster convergence to the synchronized state.

III.5 On global stability of the in-phase synchronized state

In all considered examples, the optimized interactions yield the in-phase synchronized state as a globally stable fixed point, even if the direct interaction gives bistability of the in-phase and anti-phase states (see Figs. 3, 9, and 13) footnote . Since global stability of the in-phase synchronized state is not imposed in the present optimization problem, these results are coincidental, as stated before. In the present examples, low-order Fourier components are dominant in the waveforms of the oscillation and phase sensitivity function, and thus the phase coupling function $\Gamma(\phi)$ is also a long-wave smooth function of the phase $\phi$. By maximizing $-\Gamma^{\prime}_{a}(0)$, the first pair of positive and negative lobes of $\Gamma_{a}(\phi)$ on both sides of $\phi=0$ axis is further enhanced, resulting in the global stability of the in-phase synchronized state. In the case of the nonlocal interaction, linear filtering of the field variable also contributes to the smoothness of the phase coupling function. In contrast, in the case of the nonlinear interaction, there is no such smoothing effect by the filtering and the phase coupling function is less smooth than the nonlocal case, as we can observe in Figs. 3 and 13.

IV Summary

We developed a method for optimizing the interaction function between a pair of mutually interacting reaction-diffusion systems exhibiting stable rhythmic patterns. We first considered the case with nonlocal interaction and derived the optimal filtering functions for two types of constraints. We then showed that the optimal nonlinear interaction function is given by the negative of the derivative of the phase sensitivity function in the general case. Using the FitzHugh-Nagumo reaction-diffusion system, we illustrated that the synchronization between the two systems becomes much faster than the case where all field variables are directly interacting with identical intensity. These results could be tested experimentally, e.g., by using a pair of spiral patterns in spatially extended photosensitive Belousov-Zhabotinsky reaction Hildebrand ; ref:showalter15 ; ref:weiss15 ; ref:weiss17 .

Acknowledgements.

Y.K. acknowledges financial support from JSPS (Japan) KAKENHI Grant Number JP16K17769. S.S. acknowledges financial support from JSPS (Japan) KAKENHI Grant Number JP15J12045. T.Y. acknowledges financial support from JSPS (Japan) KAKENHI Grant Number JP15K05221. H.N. acknowledges financial support from JSPS (Japan) KAKENHI Grant Numbers JP16H01538, JP16K13847, and JP17H03279.

Appendix

Each component of the optimal nonlocal interaction

[TABLE]

can be expressed as

[TABLE]

where

[TABLE]

That is, the optimal nonlocal interaction can also be represented as

[TABLE]

and we use this expression in the numerical simulations. When we discretize a $d$-dimensional spatial domain by using $N^{d}$ grid points and one period of oscillation by using $M$ points, the computational costs of Eq. (53) and Eq. (58) are $O(N^{2d})$ and $O(MN^{d})$, respectively. Therefore, Eq. (58) is much more efficient than Eq. (53) when the spatial dimension $d$ of the system is large.