# Universal Scaling of Wave Propagation Failure in Arrays of Coupled Nonlinear Cells

\[

## Abstract

We study the onset of the propagation failure of wave fronts in systems of coupled cells. We introduce a new method to analyze the scaling of the critical external field at which fronts cease to propagate, as a function of intercellular coupling. We find the universal scaling of the field throughout the range of couplings, and show that the field becomes exponentially small for large couplings. Our method is generic and applicable to a wide class of cellular dynamics in chemical, biological, and engineering systems. We confirm our results by direct numerical simulations.

\]

The impact of discreteness on the propagation of phase fronts in biophysical, chemical, and engineering systems has been intensively studied during the last decade. Among the diverse examples are calcium release waves in living cells , reaction fronts in chains of coupled chemical reactors , arrays of coupled diode resonators , and discontinuous propagation of action potential in cardiac tissue . All these disparate systems share a common phenomenon of wave front propagation failure, independently of specific details of each system. Recently this effect has drawn considerable attention (see, e.g., Refs. ). Numerous experimental evidences show that the propagation failure occurs at finite values of the coupling strength (a critical coupling. This is contrary to continuous systems, where wave fronts propagate for arbitrary couplings . A challenging problem is to establish the universal properties of the critical coupling; this is crucial for making predictions of qualitatively different regimes of system dynamics.

In this Letter we consider the universal behavior of phase separation fronts in one-dimensional nonlinear discrete systems in an external field. We study the propagation failure transition for a class of simple dynamical models describing experimental observations in arrays of coupled nonlinear cells, such as chains of bistable chemical reactors , systems of cardiac cells , etc. A new analytical method is presented to study generic properties of the critical external field of the transition. We find, using this method, how the critical field scales with the intrachain coupling. This method is applicable for a wide range of the couplings. We confirm our analytical predictions by direct numerical simulations of the full system. Our model in general is given by the following set of coupled nonlinear equations:

$$\gamma \frac{du_n}{dt}=C(u_{n+1}+u_{n1}2u_n)\frac{G(u_n,E)}{u_n}.$$
(1)   
Here $`u_n`$ is the order parameter at the $`n`$-th site, $`\gamma `$ is the damping coefficient, $`C`$ is the coupling constant, and $`G(u_n,E)`$ is the onsite potential, where $`E`$ is the applied field. The potential has at least two minima separated by a barrier $`u_B`$. The external field $`E`$ is responsible for the energy difference between the minima. This provides for one globally stable and one metastable minimum, $`u=u_+`$ and $`u=u_{}`$, respectively. Phase fronts connect the two minima and tend to propagate, to increase the size of the energetically more favorable phase, $`u_+`$.

The mechanism of front propagation is manifested in the competition between the system discreteness and the driving field $`E`$. This competition gives rise to the propagation failure at the critical field $`E_c`$, which depends upon the intercellular coupling . For $`E>E_c`$, the front propagates at the velocity $`V`$ (see, e.g., Ref. ) vanishing at the transition point, $`E=E_c`$.

The results presented in this Letter are generic and apply to a wide class of systems, independently of the details, as long as the potential $`G(u,E)`$ has the bistable structure described above . We start our analysis with the sine-Gordon potential, $`G=K(1\mathrm{cos}uEu)`$, with the factor $`K`$ being the potential amplitude. This potential was chosen in order to address the systems with phase dynamics, where the order parameter possesses a natural periodicity, such as arrays of Josephson junctions . The dimensionless dynamics in this case are given by

$$\frac{du_n}{dt}=\beta (u_{n+1}+u_{n1}2u_n)+E\mathrm{sin}u_n,$$
(2)   
where time is rescaled as $`ttK/\gamma `$ and the dimensionless coupling is $`\beta =C/K`$. Our consideration is focused on the elementary fronts connecting two nearest minima out of the infinite set of the potential minima. The system is invariant with respect to the shift $`uu+2\pi `$, so we choose without loss of generality $`u_{}=\mathrm{arcsin}E,u_+=u_{}+2\pi ,u_B=\pi \mathrm{arcsin}E`$. This makes the dynamics effectively bistable.

We introduce a novel analysis to find the universal dependence of the critical field $`E_c`$ on the coupling $`\beta `$. For not too large $`\beta `$, the results are in a good quantitative agreement with the corresponding asymptotic description. Finally we address a bistable fourth degree polynomial potential with the corresponding force in Eq. (1), $`G/u=u(u1)(u1/2+E)`$. Such a potential is applied, e.g., to describe the propagation failure in arrays of chemical reactors  and coupled diode resonators . This illustrates the applicability of our approach to generic potentials. The obtained analytical results for both potentials are justified by the numerical simulations of the full systems.

For large dimensionless coupling, $`\beta 1`$, Eqs. (2) approach the continuous regime described by the overdamped sine-Gordon equation, $`u_t=\beta u_{xx}\mathrm{sin}u+E`$. Here $`x`$ is a spatial coordinate standing for the continuous site number. This equation possesses a front solution, which propagates at nonzero velocity, for any driving field $`E`$. For small field, the front has the form $`u4\mathrm{arctan}\left[\mathrm{exp}(z/\sqrt{\beta })\right]`$, with the traveling wave coordinate $`z=xVt`$ and velocity $`VE\sqrt{\beta }`$ . Therefore, the critical field $`E_c`$ vanishes in the continuous regime. For finite $`\beta `$, however, the external driving can be balanced by the effect of discreteness, thus leading to a propagation failure at a finite value of $`E=E_c(\beta )`$. Note that for a wide class of bistable systems, there is a global instability at $`E=E_{gl}`$, above which the potential minima $`u_\pm `$ (or one of them) cease to exist. Typically the propagation failure occurs for $`E_c<E_{gl}`$. In particular, for Eqs. (2) one has $`E_{gl}=1`$.

In order to study the onset of the propagation failure, one has to consider the stationary case of Eqs. (2), i.e., $`u_n/t=0`$ (Fig 1). We denote the site closest to the barrier separating two potential minima, as the “front site” and assign to it number $`n=0`$.

Before developing the main approach of the present study, we briefly present our results obtained by means of single-active-site theory, valid for the discrete regime of not too large coupling $`\beta `$. This theory is based on the fact that, for such $`\beta `$, only the front site, $`n=0`$, experiences nonlinearity, while all the other sites are close enough to either $`u_+`$ or $`u_{}`$ to be in a linear regime, see Fig 1. Solving the linearized version of Eqs. (2) in the stationary case, one obtains $`u_n=\mathrm{\Delta }+A\mathrm{exp}(\lambda n)`$ and $`u_n=2\pi +\mathrm{\Delta }+B\mathrm{exp}(\lambda n)`$, for $`n\mathrm{\hspace{0.33em}0}`$ and $`n\mathrm{\hspace{0.33em}0}`$, respectively. Here $`\mathrm{\Delta }=\mathrm{arcsin}E`$, $`\lambda =\text{arccosh}[\mathrm{\hspace{0.17em}1}/(2\beta )+\mathrm{cos}\mathrm{\Delta }]`$. Matching these at $`n=0`$ leads to $`A=2\pi +B`$, which after substitution in the full nonlinear equation for $`u_0`$, yields

$$2\beta (\pi +B)(e^\lambda 1)+\mathrm{sin}(B+\mathrm{\Delta })=E.$$
(3)   
For $`E<E_c`$, Eq. (3) possesses two solutions, corresponding to stable and unstable fronts in the original problem (2). At the bifurcation point $`E=E_c`$ these solutions merge and disappear, so that for $`E>E_c`$ no stationary solution to Eqs. (2) exists, and the front propagates. This implies that at the bifurcation point the derivative of the left hand side of Eq. (3) with respect to $`B`$ must equal zero. Then we obtain, after some calculations, an approximate expression for $`E_c`$:

$`E_c`$ $`{\displaystyle \frac{\sqrt{2\sqrt{1+4\beta }4\beta 1}f(\beta )\mathrm{arccos}[f(\beta )]}{1+f(\beta )}},`$ (5)   
$`\text{where}f(\beta )=2\left[1{\displaystyle \frac{\beta +1}{1+2\beta +\sqrt{1+4\beta }}}\right].`$   
The graph of $`E_c(\beta )`$ given by (5) is shown in Fig. 2 (dotted line). We see from the figure that it agrees well with the results of numerical simulations of the full system (2), for small to moderate values of $`\beta `$ ($`\beta 0.5`$). This makes the single-active-site theory substantially more useful than a regular small $`\beta `$ perturbation theory, which works only for much smaller $`\beta `$ ($`0.1`$).

Although the single-active-site theory gives reasonable predictions for the propagation failure transition for moderate values of the intrachain coupling $`\beta `$, it does not provide a universal scaling of the critical field $`E_c`$ with $`\beta `$ (see Fig. 2). This motivates developing a general theory to describe the phenomenon of the propagation failure throughout the range of $`\beta `$. Such a theory is presented below.

A stationary front in Eq. (2) can be obtained as an extremum of the free energy

$$=\underset{n}{}[\frac{\beta }{2}(u_{n+1}u_n)^2+1\mathrm{cos}u_nEu_n].$$
(6)   
For the case of zero field $`E`$, the front is always stationary, and no transition exists. In this case we have for the front, taking into account the effect of discreteness (see also )

$`u_n=w(n){\displaystyle \frac{\mathrm{sin}w(n)}{12\beta }},`$ (7)   
$`w(x)=2\mathrm{arctan}\left[b\mathrm{sinh}^1\left({\displaystyle \frac{bx}{\sqrt{\beta }}}\right)\right],x0,`$ (8)   
$`w(x)=2\pi w(x),x0,`$ (9)   
with the coefficient $`b=[\mathrm{\hspace{0.17em}1}1/(12\beta )]^{1/2}`$.

At zero field $`E`$, the free energy $``$ (6) possesses an infinite set of minima of equal depths, separated by barriers. Each of these minima corresponds to a stable stationary front and differs from the other fronts by a shift on an integer number of sites. Each barrier corresponds to an unstable front. When the external field $`E`$ is applied, the set of minima is tilted, so that the difference between their depths is determined by the field value. As the field increases, each of the minima approaches its adjacent barrier (down the energy landscape), until they finally merge, at $`E=E_c`$. For $`E>E_c`$, the energy has no extrema left, and therefore the fronts propagate without limit.

Our method is based on choosing profiles $`\{u_n\}`$ that can be parameterized by an effective coordinate $`\alpha `$, and replacing the argument $`n`$ of function $`w(n)`$ in (7) with $`n+\alpha `$. We show below that this choice allows us to find the critical field $`E_c`$. It can be demonstrated that, for $`E=0`$, the points $`\alpha =1/2+m`$ ($`m`$ is an integer) correspond to stable fronts, and $`\alpha =m`$, to unstable fronts. The former and the latter are energy minima and maxima, respectively, satisfying $`d/d\alpha =0`$. Note that the chosen parameterization by $`\alpha `$ may be interpreted as a continuous shift of a profile $`\{u_n\}`$.

Our goal now is to evaluate the free energy landscape $`(\alpha )`$ given by (6), for nonzero field $`E`$. When $`E=E_c`$ the minima of $`(\alpha )`$ merge with their adjacent maxima at the inflection points $`d^2/d\alpha ^2=0`$. We have substituted the fronts (7) with $`w(n+\alpha )`$ into energy (6) and found $`E_c`$ for various $`\beta `$’s, using Mathematica software. The result is shown in Fig. 2 (dashed line).

In order to find the analytical form of $`(\alpha )`$, we use the Poisson summation formula for infinite series:

$$\underset{\mathrm{}}{\overset{\mathrm{}}{}}F(n)=_{\mathrm{}}^{\mathrm{}}F(x)\left[\mathrm{\hspace{0.17em}1}+2\underset{k=1}{\overset{\mathrm{}}{}}\mathrm{cos}(2\pi kx)\right]𝑑x,$$
(10)   
where $`F(n)`$ is the $`n`$th element of the series in (6). One finds that the integral $`F(x)𝑑xE\alpha `$, up to a term independent of $`\alpha `$. Then, after closing the integration contour in the complex plane, we see that $`k>1`$ terms in (10) are exponentially small, compared to the $`k=1`$ term. This dominant term can be evaluated approximately to give the following expression for $`(\alpha )`$

$$(\alpha )=2\pi E\alpha +\mathrm{\Omega }\beta \mathrm{cos}(2\pi \alpha )\mathrm{exp}\left[\frac{\pi ^2\beta }{\sqrt{\beta 1/12}}\right].$$
(11)   
The factor $`\mathrm{\Omega }`$ is a constant and can be found by comparison of (11) with the results obtained by Mathematica (see Fig. 2), which yields $`\mathrm{\Omega }340`$. Then, using the conditions $`d/d\alpha =d^2/d\alpha ^2=0`$, we find

$$\alpha =\frac{1}{4},E_c=\mathrm{\Omega }\beta \mathrm{exp}\left[\frac{\pi ^2\beta }{\sqrt{\beta 1/12}}\right].$$
(12)   
The general scaling $`E_c(\beta )`$ of the propagation failure transition given by (12), virtually coincides with the one obtained by Mathematica, shown in Fig. 2 . In the continuous regime of large $`\beta `$ the result (12) has the following simple asymptotic form $`E_c=\mathrm{\Omega }\beta \mathrm{exp}[\pi ^2\sqrt{\beta }]`$. This implies that the critical field $`E_c`$ decays exponentially for large couplings $`\beta `$.

To confirm our general theory, we have performed numerical simulations of the full system (2). We have used the implicit second order integration method, in order to obtain the transition line $`E_c(\beta )`$. The results are presented in Fig. 2 (solid line). We see in the figure a good quantitative agreement of our analytical predictions with the results of the simulations.

To address the propagation failure in arrays of chemical reactors  and diode resonators , we turn to a fourth degree polynomial potential leading to a cubic force in Eqs. (1): $`G/u=u(u1)(u1/2+E)`$. The globally stable and metastable minima for this potential are $`u_+=1`$ and $`u_{}=0`$, respectively; the barrier $`u_B=1/2E`$. Applying our general theory developed above, we find the propagation failure transition line $`E_c(\beta )`$. In Fig. 3 we show the results obtained by Mathematica (dashed line). We see again that the theory is confirmed by numerical simulations of the full system (solid line). The analytic expression of $`E_c(\beta )`$ analogous to (12) is cumbersome , so in this Letter we only give its asymptotic form for large $`\beta `$: $`E_c=\mathrm{\Omega }\beta \mathrm{exp}(\eta \sqrt{\beta })`$. This form has the same structure as the asymptotics of critical curve (12) for sine-Gordon potential, which provides a persuasive argument that the proposed approach is generic, independent of the specific potential. The values of $`\mathrm{\Omega }`$ and $`\eta `$ are found to be $`\mathrm{\Omega }429,\eta 21.9`$.

To conclude, we have introduced a general method to analyze the transition of front propagation failure in arrays of coupled nonlinear cells. Using this method we have determined, for bistable dynamics of the cells, how the critical external field of the transition scales with intrachain coupling. To demonstrate the generality of the new method, we have carried out our analysis for two different effectively bistable potentials: (i) the sine-Gordon potential, useful e.g., for describing arrays of Josephson junctions ; and (ii) a fourth degree polynomial potential applicable to chains of chemical reactors  and coupled diode resonators . Our theoretical predictions have been confirmed by numerical simulations of the full systems. After this study of the propagation failure problem in the one-component systems (1), the following question arises: what is the mechanism of this phenomenon in systems described by more complicated dynamics? For example, the “single-pool-model” for intracellular calcium waves  is represented by two dynamic components: the calcium concentration and the fraction of vacant ionic channels. An important and virtually unexplored issue is how the effect of discreteness of ionic channels determines the onset of the propagation failure of calcium waves. We expect that our techniques can provide new insights into this and other types of more complicated chains of coupled nonlinear elements in biophysics, chemistry, and engineering. A further challenging direction of study here is to extend our analysis of the propagation failure to two- and three-dimensional systems. This would allow one to understand the effect of discreteness on the dynamics and stability of such structures as spiral and scroll waves in cardiac tissue .

We thank John Pearson for fruitful discussions and valuable advice. This work was supported by the Department of Energy under contract W-7405-ENG-36.