Lateral inhibition in relaxation oscillators provides a basis for computation

Complex spatio-temporal patterns seen in natural phenomena, ranging from the intricate symmetric geometry of snowflakes [1] to the development of form in multicellular organisms [2, 3, 4] and large-scale plankton patchiness in oceanic algal blooms [5, 6, 7], are associated with processes that are typically far from equilibrium [8, 9]. The behavior of the dissipative systems where such patterns arise can be described using a dynamical framework, yielding solutions ranging from waves, breathers, and solitons, to chaos [10, 11]. Each of the collective states, which are characterized by distinct emergent phenomena, encodes information about the instantaneous configuration of the system. Consequently, a transition from one pattern to another can be viewed as a “computation” that results in the transformation of such information [12, 13, 14]. Indeed, the dynamical evolution of systems such as spin glasses, which flow down rugged energy landscapes to one of many possible local minima [15], have been used to implement neural network-inspired computation [16, 17]. This paradigm has been successfully used to model brain function [18] as well as to provide innovative solutions to hard combinatorial optimization problems [19, 20, 21]. However, it has been challenging to formulate such a language, which connects dynamics and computation, for describing systems that are far from equilibrium.

In this paper, we present a framework for uncovering the computational logic of systems undergoing transitions between multiple non-equilibrium steady states. Two archetypes of the irreversible processes that characterize far-from-equilibrium systems are diffusion and chemical reactions (such as combustion) occurring in open systems [22]. Apart from being fundamental natural processes, they play important roles in biology, for instance, allowing transportation of molecules via diffusion, and their synthesis and replication through metabolic reactions [23]. Reaction-diffusion systems, in which spatio-temporal patterns can emerge through self-organization [24, 25], provide a natural framework for describing information processing in terms of out-of-equilibrium dynamics. Indeed, systems embodying these principles, such as oscillating chemical reactions and slime molds, have been used to implement specific computational algorithms [26, 27, 28, 29, 30, 31, 32]. Computation in these systems is implemented by the dynamical evolution of the initial (input) state to the final (output) state. The general computational logic presented here suggests that these algorithms belong to a much broader paradigm, in which such evolutions can be completely specified by a set of transformation rules (logic) that map the set of non-equilibrium steady states of the system to itself. By altering the perturbations applied to the system, the logic can be changed, allowing different types of computations to be realized (analogous to programming a computer). This can potentially address a central problem in biology: how living organisms “compute”, i.e., process information vital for survival through biochemical means [33, 34].

As chemical reactions occurring inside living organisms exhibit temporal recurrence across scales, we illustrate this paradigm using a generic non-equilibrium system of relaxation oscillators that are coupled diffusively. The resulting collective dynamics generates a variety of spatiotemporal patterns depending on kinetic rates and coupling strengths [35, 36, 37]. Unlike previous attempts at implementing computation in a reaction-diffusion framework, which typically involve interactions between traveling waves [38, 33, 39, 40], we use a temporally invariant, spatially heterogeneous non-equilibrium steady state. Such Spatially Patterned Oscillator Death (SPOD) states are observed for relatively strong coupling between oscillators [36]. The temporal invariance allows a particularly succinct dynamical description, and we show that a potentially infinite number of phase space trajectories resulting from different initial conditions reduce to a finite number of qualitatively distinct states. This allows coarse-graining of the state space continuum to a discrete set that can be mapped to distinct binary sequences. The time-evolution of the system thereby reduces to operations in which symbolic strings are transformed to one another. In analogy with computation in a Turing machine [41], here binary strings representing the initial and final states are the input and output, respectively, while the role of the program is played by the transformation-inducing perturbation, viz., the binary pattern of stimuli (i.e., on/off) applied to the array of oscillators. We generate a comprehensive taxonomy of the states, investigate their local, as well as, global stabilities, and provide a complete catalogue of perturbation-induced transitions between them, yielding a fundamental template for realizing computation through the collective dynamics of non-equilibrium systems.

To describe the dynamics of each node in the system of coupled relaxation oscillators, we have used a generic model for such oscillators, viz., the FitzHugh-Nagumo (FHN) equations. This model comprises a fast activation variable $u$ and a slow inactivation variable $v$ evolving as

where the parameters $\alpha(=0.139)$ and $k(=0.6)$ describe the kinetics, $\epsilon(=0.001)$ governs the recovery rate and $b(=0.16)$ characterizes the asymmetry of the limit cycle dynamics. The simplest coupling topology that allows the system to be used for computation is that of a one-dimensional ($1$-D) array, as each distinct state can be mapped to a symbolic sequence by using appropriate thresholds. In this paper we have investigated a $1$-D system of $N$ coupled oscillators with different boundary conditions, corresponding to a chain and a ring. Following Ref. [36], we consider the oscillators to be coupled by diffusion exclusively via the inactivation variable $v$. Such a coupling captures the essence of the physical process underlying experiments in microfluidic devices [35], where the beads containing the oscillating chemical reactants are suspended in an inert medium that allows only the inhibitory species to diffuse. A linear chain with passive elements at the boundary [Fig. 1 (a)] approximates the geometry of such experimental systems, taking into account the presence of the inert medium [36]. The dynamics of the resultant system is described by:

where $i=1,2,...,N$ ($N=10$, unless specified otherwise). The boundary conditions are specified by $\dot{v}_{0}=D_{v}(v_{1}-v_{0})\,,$ $\dot{v}_{N+1}=D_{v}(v_{N}-v_{N+1})\,.$ The diffusion coefficient $D_{v}(=2\times 10^{-3})$ denotes the strength of coupling between the inactivation variables of neighboring relaxation oscillators [42].

The effect of external stimulation on a node $i$ of the system, such as optical illumination in the case of light-sensitive chemical reactions that can initiate oscillatory activity in the medium [43], is incorporated through the term $A_{i}$ which is assumed to be constant for the duration of stimulation. In the simulations reported here, the signal intensity $A_{i}$ is assumed to be the same ($=A$) for all nodes. Varying $A$ and the stimulation duration $\tau$ allows us to investigate the effect of signal strength on the collective dynamics. A stimulation protocol is defined by the choice of nodes $i$ (=$1,\ldots,N$) for which $A_{i}\neq 0$ corresponding to a total of $2^{N}-1$ possible protocols. In principle, these protocols can implement logic gates. Indeed, we explicitly show that a conservative logic gate [44] that is capable of universal computation, can be realized using such protocols [45]. We have verified that the results reported in this paper are robust with respect to different choices of $N$, as well as, small changes in $b,\alpha,k,\epsilon$, $D_{v}$ and $A$.

Fig. 1 (b) shows the system (2) converging to a characteristic SPOD state, marked by a temporally invariant pattern of high and low values for $u_{i}$, as well as, $v_{i}$ ($i=1,\ldots,N$). Note that, a low (high) value of $u$ implies a low (high) value of $v$ for each element, and vice versa. For each SPOD state, all elements that are low have exactly the same numerical value for $u$ (and $v$), as do all elements that are high. The configuration can hence be mapped to a binary sequence using an appropriate threshold. As the steady-state probability distributions of $u,v$ (estimated from $10^{4}$ realizations) have a clearly bimodal form, we choose a threshold (viz., $u=0.4,v=0.06$) from the region between the peaks, our results being robust with respect to this choice. To allow an arbitrary binary string to be used as the initial state of the system [e.g., as in Fig. 1 (c)], we map the entries ($0,1$) of the string to the $u,v$ values that occur at the two peaks of the distributions (viz., $0:u=-0.1,v=0.03$ and $1:u=0.85,v=0.093$).

The system exhibits multistability, i.e., for a given choice of parameter values, it can converge to any one of multiple possible SPOD states, each having a corresponding binary representation. Fig. 1 (c) shows an example of how one such binary string can be transformed to another through an external perturbation - a property that is vital for the system to be used for computation. While, in principle, a system of $N$ elements should be able to represent $2^{N}$ binary strings, only a small fraction of these can actually be obtained. For example, we observe only $68$ distinct strings for a chain of $N=10$ oscillators, and $122$ distinct strings for a ring of the same size (which reduces to $14$ on considering the cyclic permutations to be identical, taking into account the periodic boundary condition). This restriction arises from the physical mechanism by which SPOD states are generated, and can be understood by considering a pair of coupled oscillators. When the two oscillators are in opposite branches [Fig. 1 (d), right], the forces arising from diffusion and intrinsic kinetics acting on the oscillators are balanced under strong coupling, resulting in either a $0$-$1$ or a $1$-$0$ state [36]. For a $1$-D array of $N$ coupled oscillators, we observe that the SPOD states contain at most two consecutive 0s or 1s, hence limiting the number of observable binary strings.

The instability of sequences containing long ($>2$) contiguous blocks of 0s or 1s, for the range of $D_{v}$ we use, arises because the diffusion term for the interior node(s) in such a block is negligible (thereby disrupting the balance of forces responsible for the arrest of oscillations). Thus, these nodes effectively become uncoupled and hence oscillate, destabilizing any time-invariant states containing such blocks [45]. This also implies that SPOD states that begin or end with ‘00’ or ‘11’ will not be observed in a chain with passive elements at its ends [45]. We note also that the inversion of an observed SPOD state, obtained by exchanging $0$s with $1$s and vice versa, yields another allowed SPOD state. The stability of two coupled oscillators arrested in a low and a high activity state suggests that if such pairs can completely tile the entire sequence represented by the array [Fig. 1 (d), left], the pattern will be stable. Indeed, we observe that all SPOD states observed in system (2) can be viewed as a sequence of stable pairs of adjacent oscillators, each pair comprising one oscillator arrested in a high and the other in a low activity state (i.e., ‘01’ or ‘10’). This provides a basis for viewing the system as an anti-ferromagnetic spin chain (spin up/down corresponding to high/low activity), allowing us to formulate an effective energy function for the system (see below).

A simple combinatorial argument allows us to obtain the number of distinct SPOD states appearing in a chain of $N$ oscillators, $n_{SPOD}(N)=2F_{N}$, the $N$-th term of the Fibonacci sequence [45]. For a ring, a more involved combinatorial argument is required to determine the distinct number of SPOD states, which in this case is the number of cyclic strings with at most two consecutive 0s and 1s anywhere in the string. For a system of size $N$, this number is the coefficient $C(N)$ of the generating function $f(s)=\sum C(N)s^{N}=1+2s+4s^{2}+6s^{3}+6s^{4}+O(s^{5}),$ obtained using the Goulden-Jackson cluster method [46, 45]. Thus, $C(0)=1$, $C(1)=2$, $C(2)=4$ and for $N>2$, $C(N)=\phi_{+}^{N}+\phi_{-}^{N}+\omega^{-N}+\omega^{N},$ where $\phi_{\pm}=\frac{1\pm\sqrt{5}}{2}$ and $\omega=\sqrt[3]{1}$. Taking into account the periodic boundary conditions of the ring and using Polya enumeration theory [46], the number of strings unique under cyclic permutations is obtained as $UC(N)=\frac{1}{N}\sum_{d:d|N}\varphi(d)C^{\prime}(\frac{N}{d})\,,$ where $C^{\prime}(1)=0,C^{\prime}(2)=2,$ and $C^{\prime}(j)=C(j),$ for $j>2$. Here $\varphi(n)$ is the Euler’s totient function (or phi function) defined as the number of positive integers $\leq~{}n$ and relatively prime to $n$. A numerical evaluation of these formulae show that the number of distinct SPOD states increases exponentially with system size $N$, independent of the connection topology [45].

We now investigate the effect of perturbing the SPOD states that occur in a ring of $N$ elements by using each of the $2^{N}-1$ possible stimulation protocols. Fig. 2 (a) shows the perturbation-induced transitions between strings corresponding to the $14$ distinct states (invariant under cyclic permutations) for $N=10$. The number of $0$s present in a string indicates its relative stability, measured in terms of the number of transitions that originate from, or terminate at, the string. In particular, strings with the maximum (minimum) allowed number of $0$s, i.e., $6$ ($4$) for $N=10$, act as sinks (sources, respectively), although sources and sinks having $5$ zeros do also occur in this case. The maximum number $n_{\max}(0)$ of $0$s allowed in a binary string of length $N$, subject to the constraint that substrings of contiguous 0s (or 1s) of length $>2$ are not allowed, is equal to the largest integer $\leq N/3$, while the minimum number $n_{\min}(0)$ of $0$s allowed is $N-n_{\max}(0)$. Note that the strings that are most frequently obtained when starting from random initial conditions [the largest nodes in Fig. 2 (a)] are not necessarily the ones most stable to perturbation. This difference between (i) the probability $P_{\rm pert}$ that a particular string results from perturbing the allowed states, and (ii) its basin size, i.e., the probability of obtaining the string from random initial conditions, can be understood in terms of the distinction between local and global stabilities.

The local stability of a string can be quantified in terms of the decay rate of perturbations to the corresponding SPOD state, which is related to the largest eigenvalue $\lambda_{\max}$ of the Jacobian obtained by linearizing Eqn. (2) around this state. As seen from Fig. 2 (b), strings that have the highest $P_{\rm pert}$ (and hence are sinks), are associated with the most negative values of $\lambda_{\max}$, implying that they have much higher local stability compared to sources, i.e., strings having the lowest $P_{\rm pert}$. The global stability of each string can be defined in terms of an energy function using the mapping to an anti-ferromagnetic spin chain. Identifying the $1$ ($0$) in the $i$th node of a string with $S_{i}=+1$ or “up” ($S_{i}=-1$ or “down”, respectively) for the corresponding spin, the associated energy is defined as $E=\sum_{i=1}^{N}S_{i}S_{i-1}$ (periodic boundaries imply $S_{0}=S_{N}$). Thus, states containing adjacent spins in the same orientation would have a higher energy, decreasing the likelihood that they will emerge from an arbitrary initial state and hence reducing their basin size [Fig. 2 (c)]. The use of an energy landscape to describe the system dynamics is consistent with trajectories initiated from each of the states converging to absorbing (sink) states, whose existence implies the absence of periodic cycles.

To conclude, the results presented here point to a new computing paradigm involving the global nonlinear response of a continuous medium to local perturbations, which is coordinated through diffusive transport. Our results could be verified in an experimental set-up with photoperturbations applied to arrays of light-sensitive chemical droplets in microfluidic devices [35, 43, 47]. As our model is generic, it may also be possible to execute such a computation scheme in other experimental systems involving coupled relaxation oscillators [48], including engineered bioelectric tissues [49], phototactic Physarum plasmodium [50] and insulator-metal state transition devices [51]. Considering strong inhibition between neighboring units allows us to use time-invariant patterns as the fundamental states in computation involving reaction-diffusion, instead of relying on interactions between traveling waves [40]. As there are many biological processes which involve lateral inhibition between elements having rhythmic activity that can be modeled as relaxation oscillators [52], the framework we present here may offer insights into how computation may be implemented in living systems [53].