Modeling somatic computation with non-neural bioelectric networks

(b) The overall performance of the best evolved pattern detector, with data collected from a set of 1000 simulations: 100 parallel sets each with a random sequence of 10 simulations. As expected, the pattern detector classifies patterns similar to the French flag (a total of 500 sample inputs) as “French flag”, and those that are dissimilar (500 inputs) as “not French flag”. This suggests that the pattern detector is robust to noise to a sensible extent. The width of the classification boundary was set implicitly due to the way sample input patterns from the two classes were generated (details in the ‘Methods’ section). (c) The behavior of the best French-flag detector shown for a set of four representative cells (three inputs out of a total of eighteen and one output) for a random sequence of five patterns. Inset highlights the four nodes whose colors correspond with those in the time series. Highlighted in grey squares are the cases where a French-flag-like pattern is input for which the output is depolarized, as expected; for all other cases where a random pattern is shown, the output is hyperpolarized.

We conclude that BEN has in principle the ability to distinguish between patterns. The biological implication is that somatic tissues can in principle have the same ability, and thus be a part of much larger pattern regulation mechanisms.

Biological processes are complex by nature. The rampant complexity is partly managed by nature by way of modularity78, in gene regulation79 and the brain80,81, for example, at multiple scales. Modularity is beneficial because the modules can be independently tinkered with for the purpose of large-scale outcomes79. Thus, it is important to understand whether bioelectric circuits can implement logic functions that are too complex to be described by a single gate but are actually combinations of elementary logic primitives.

Compound logic gates are compositions of modules of logic gates, and their biological equivalents underlie several processes including sporulation in B. subtilis and the neuronal dynamics of C. elegans82,83. They can also have pharmaceutical applications—for example, the control of bacterial invasion of tumor cells and diagnosis of cellular environments and automatic release of drugs82.

One way to design compound networks, in general, is to compose them from appropriate pre-designed modules. For example, a “NAND” gate can be constructed by combining the AND and NOT gates. Compound genetic circuits have been constructed in this fashion84,85, as well as designed de novo from scratch82. Here, we show how to construct compound logic gates in physiological networks by composing pre-trained BEN logic gate modules. Unlike genetic networks and neural networks, which are directed, BEN networks are bidirectional, due to the symmetric nature of gap junctions. Thus, unlike genetic networks, compound logic gates in BEN cannot be constructed simply by connecting the modules. Put more precisely, if the modules were simply connected together, then the output of one module will also receive signals from the input of the connected module. This could potentially result in the upstream module outputting the wrong state, and thus the downstream module receiving and outputting the wrong states as well. To mitigate this, we included a “bridge” layer that interfaces the modules and trained it so that the compound gate as a whole behaves as desired: the downstream module outputs the correct state for every set of inputs received by the upstream module. For more details on the architecture and training of the bridge, see Supplementary 4. Figure 7 shows an example of a successfully trained NAND gate by composing pre-trained AND and NOT gates.Figure 7The behavior of a successfully trained NAND gate for a random sequence of four inputs. Orange and purple lines represent the activities of the inputs, while green represents the output. The vertical dashed lines mark the beginning of each trial in the sequence, however the inputs are applied, after an initial transient, at the time points marked by the grey triangles.

Inset is a schematic of the architecture of the gate: pretrained AND and NOT gates are connected by a “bridge”, a single layer of cells, that is alone trained.

We conclude that it is possible to construct compound logic gates by composing pre-trained modules and training only an interface bridge connecting the modules. This method can in principle be extended to more complex circuits.

Various forms of natural and artificial computing systems with unconventional mechanisms have been proposed86. Even the modern von Neumann computer architecture was inspired by the mechanisms of DNA transcription and translation87. Other examples include “liquid computers” that involve interactions of fluid streams and droplet flows to perform logic operations88; “chemical computers” that utilize the dynamics and equilibrium-behavior of chemical reactions89–93 for logic operations and learning94,95; DNA-based computers that utilize various forms of catalytic reactions such as ‘seesawing’ and ‘strand-displacement’ for pattern recognition96, protein-DNA interactions to implement logic functions97 and even for general-purpose computing98; gene regulatory networks that take advantage of the multistability of gene activity for learning99, memory100–102 and logic84,85,101,102; “bacterial computers” that use enzyme-based computing to implement logic gates103; electric current flow making use of Kirchhoff’s law for solving mazes104; and even “sandstone computers” that utilize natural erosion processes to implement logic gates105. BEN is an unconventional computing system, since it is a non-neural system like the examples above.

Our work does not require a commitment to the view that every biological question must be dealt with via a computational perspective, despite the fact that there is a robust body of new experimental work driven by the hypothesis that some tissues are indeed implementing processes we would recognize as memory, comparison, or measurement53,106–108. The key question we answer here is: are pre-neural physiological networks able, in principle, to perform simple computations such as logical operations? Prior to these results, it was not known, and many developmental biologists’ intuitions are that such networks cannot perform such functions - they are almost never considered for model-building. Our data provide the first proof-of-existence showing that the class of biological systems where a computational approach could be useful, such as evolved or synthetic regulative morphogenesis, needs to be expanded to include somatic bioelectric networks.

BEN is also a type of patterning system. Patterning is the ability of a system to generate, sustain, and recognize patterns in the states, whether physical, chemical, or biological, of the system. It is a widespread phenomenon reported in biological systems, including pattern recognition of pathogen molecules by cellular receptors109, axial patterning in planaria110, spatial patterning during neural patterning111, cell polarity112, genetic patterning in bacterial populations113 and voltage patterning in neural tubes49. Here we focus on bioelectric non-neural networks of patterning. Various non-neural network models of patterning have been proposed. Some examples include the “packet-routing” system that uses an internet-like intercellular messaging system114, auto-associative dynamical systems115, and cellular automata116 based systems that combine local and global communication to model pattern regeneration. Our laboratory has produced sophisticated non-neural bioelectric network models of tissue patterning, like the Bioelectric Tissue Simulation Engine (BETSE) that incorporates detailed biophysical mechanisms of ion channels and transport processes47,51,117, as well as relatively simpler models that incorporate only the high-level functionality of the bioelectric components48,118,119. These models have been used to show that networks of ordinary tissues can give rise to gradient-patterns47,51, Turing patterns21,119, as well as oscillatory patterns118,120 that are thought to be essential for long-distance communication.

BEN differs from these models on two fronts: (1) it is biophysically simpler than BETSE, serving as a minimal model of dynamics sufficient for computation; and (2) it is more realistic than the model of Reference118, which uses “equivalent circuits” where the equilibrium Vmem is explicitly specified in the equations, whereas in BEN the equilibrium Vmem emerges from the ion channel parameters.

We have shown here that BEN networks can implement elementary logic gates, more complex tissue-level logic gates, pattern detectors and compound logic gates. Furthermore, we have demonstrated that logic can be implemented in circuits with bidirectional connections that is typical of non-neural tissues, in contrast to the conventional directed circuits like neural networks and digital electronic circuits. This implies that even though non-neural tissues may find it harder to implement formal logic (due to lower control over information flow; see Supplementary Fig. S2b, for example, where information flow is recurrent even in a small network, making it hard to control), they can achieve it. Not only can BEN networks compute, but they can also be robust to damage. One of the most important, and heretofore unexplained, properties of biological control circuits is that they continue to function computationally under significant perturbation (as seen by work on metamorphosis, stability of memory and behavioral programs during drastic tissue remodeling, and invariance of morphogenesis to changes in cell number76,121,122). BEN networks can exhibit a similar behavior, retaining function even after the removal of cells. We found that BEN achieves this by distributing information in the network (Supplementary 2).

The particular methods by which BEN networks were trained is not of central importance in this work and do not impact the main result – that BENs can be readily found which implement logic functions. We used BPTT simply as a way to discover the parameters of (train) BEN networks that perform specific logic functions or recognize patterns. In the biological world, training (learning) may happen in other ways that may be broadly classified as offline and online. Offline learning happens outside of the lifetime of an individual and at a population level, while online learning happens during the lifetime of an individual. Evolution is an example of the former, while Hebbian learning is an example of the latter. BPTT may be classified as a form of online learning, but it has been claimed to be not biologically realistic123,124. More biologically realistic forms of online learning have been proposed, for example, forward propagation of eligibility traces123, extended Hebbian learning125, dynamic self-organizing maps126 and instantaneously trained neural networks127. Such methods could potentially advance biologically realistic training of BEN networks, which we will develop in future research, but they are not necessary to find examples demonstrating that BEN networks can compute.

BEN models the bioelectric principles of a generic bioelectric network, and not all the physiological details, as these can vary widely across tissues. Individual cell types can be modeled in future effort simply by altering the ion channel parameters while retaining the basic form of the equations. Although BEN is a model of a non-neural bioelectric network, it has certain features that resemble those of a neural network, as described above (the weight and bias parameters, and the two-level nonlinear signal transformation). This may suggest that neural-like features are necessary for computation even in a non-neural network. On the other hand, it also suggests that non-neural features may support neural computation. Indeed, this is the subject matter of ongoing debates in neuroscience about how much glia contribute to neural computation128,129.

Our work contributes to the field of basal cognition by describing how a simple network of aneural cells may perform both basic and complex logic operations. However, our goal is not to validate basal cognition as a field – it is just one context (besides others, such as developmental neurobiology, evolutionary cognitive science, and synthetic bioengineering) which our findings impact. Even if some computational views of specific biological systems turn out to be wrong, our results provide value in showing how physiological circuits can be made with specific (and useful) behaviors, both for synthetic biology applications and to understand functional aspects of developmental physiology, the evolution of brain circuits, etc. In prior work, a computational perspective of cell behavior14,16,130 allowed our lab and other groups to formulate (1) ways of experimentally controlling large-scale anatomical outcomes13, and (2) formulate and solve the “inverse problem”14 of how the cells might actually achieve specific anatomical configurations, given cellular constraints and a range of starting conditions76,131. Our current work addresses (2), by demonstrating a possible approach that cells could use to implement logic functions under realistic cellular constraints.

Overall, BEN has the potential to shed light on as-yet-unexplained non-neural bioelectric information-processing. For example, when certain planarian species (highly regenerative animals) are cut they sometimes develop two heads, instead of a head and a tail, as a response to a temporary exposure to 1-octanol, a gap junction blocker. Without the gap junction blocker, they always develop a wild type morphology consisting of a single head and a tail. It has been shown that this decision to switch morphologies is bioelectrically controlled and non-neural in nature24,132. We hypothesize that the planarian bioelectric system encodes a pattern that either consists of a record of the gap junction blocker or not, and then maps it into different morphological decisions (wild-type or double-head) in a systematic way. We currently have a prototype BEN model that reproduces this phenomenon in a qualitative way, which we plan to report in the future. Finally, we have demonstrated in this work that the slower (reason described in ‘Methods’), continuous mode of non-neural bioelectric signaling, compared to the faster, pulsating mode of neural signaling, is sufficient to perform logical computations usually associated with the brain. This firmly supports the possibility of somatic computation in real biological systems. How it compares to the timescales of genetic signaling and how it could facilitate bioelectricity to be upstream of genetics53,103 is an open question that our future investigations will answer.

Our work paves the foundation for future research in training actual biological tissues for somatic computation. Ongoing research in our laboratory has identified ways by which gap junctions can be manipulated and controlled, further opening up the possibilities on this front. Overall, our research provides the conceptual and modeling foundations to understand and manipulate development and regeneration, and to construct computational synthetic living machines.

The equations that define a BEN model are shown in Fig. 8 using an example 2-cell system. The equations are shown for cell-2 only; the same equations with appropriate change in subscripts apply to cell-1 and likewise to any other cell in larger networks. The processes that the equations represent, like electrodiffusion, reaction-diffusion and gating, are marked appropriately in the figure, thus self-explanatory. The following remarks serve to better illustrate the model. As described in the introductory section, electrodiffusion utilizes two types of gradients, namely voltage and concentration, to drive the flux. Those gradients are represented by the terms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v}_{1}(t)-{v}_{2}(t)$$\end{document}v1(t)−v2(t) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${c}_{1}(t)-{c}_{2}(t)$$\end{document}c1(t)−c2(t) respectively in the electrodiffusion equations. The ion-pump equations are approximate versions of the full mechanism modeled in equations (23–29) of BETSE51. BEN makes use of just the terms involving the ion concentrations in equation (27) of BETSE51, as it lacks the other biochemical components like ATP, and by choosing appropriate term coefficients to compensate for the approximation it effectively preserves the relationship between ion concentrations and the pump flux.Figure 8Mathematical details of the BEN model. (Top) The equations that describe the dynamics in a 2-cell system illustrated in the figure on the left.

(Bottom) A glossary of variables used in the equations, the constants and initial conditions.

The dynamics of the signaling molecule is viewed as a nonlinear reaction-diffusion process as it can be abstractly represented as: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R(D(\nabla c))$$\end{document}R(D(∇c)), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla c$$\end{document}∇c is the concentration gradient, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D(.)$$\end{document}D(.) represents the 1st layer of sigmoid approximating a nonlinear diffusion, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R(.)$$\end{document}R(.) represents the 2nd layer of sigmoid approximating a nonlinear reaction.

Moreover, the nested form of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R(D(.))$$\end{document}R(D(.)) departs from the conventional form of reaction-diffusion which is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dc}{dt}=D(c)+R(c)$$\end{document}dcdt=D(c)+R(c), for example, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dc}{dt}=\nabla c+{c}^{2}$$\end{document}dcdt=∇c+c2. Nevertheless, it qualifies as reaction-diffusion since it clearly involves diffusion (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nabla c$$\end{document}∇c)and a nonlinear component (the sigmoid) that alters the net mass of c, representing production or decay depending on the sign of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dc}{dt}$$\end{document}dcdt.

Another important point to note is that even though sigmoid functions, which define the signal flux, also constitute the activation functions of neurons in neural network models38,64, there is a crucial difference: while neuron models are input-activated since neuronal synaptic junctions are directional, a cell in BEN is gradient-activated since gap-junctions are bidirectional. This is also partly the reason why non-neuronal cells tend to be slower than neural cells—diffusion is a slower process than directed currents.

There’s no set range of Vmem levels that a given BEN network will not exceed (as there are no explicit bounds on the ion concentration levels), but there are typically set to lie within [−80 mV, +80 mV]; see for example Figs. 2d, 3d, 4, 6 and 7. The concentration of the signaling molecule, on the other hand, is forced to lie with [0, 1], that is, within the generic minimum and maximum possible values of the actual concentration which is not explicitly considered in BEN.

The software to simulate, train and evolve BEN networks is available at: https://gitlab.com/smanicka/BEN.

We used the standard Euler integration method for integration the equations. We used a time step of 0.01 for the bioelectric dynamics and 0.02 for the signaling molecule dynamics. Each simulation was set up to run for a number of steps ranging from 300 to 600; smaller networks were run for fewer steps. Thus, the total simulation ranged between 300 * 0.01 = 3 simulation seconds to 600 * 0.01 = 6 simulation seconds.

Regardless of the gate, a training run was set up as follows. A single BEN network was instantiated with randomized parameters. The initial weight parameters were chosen in the range [−1, 1], and the biases in the range [0, 1]. It was then fed with a batch of four inputs and was simulated for about 300 time-steps, with the inputs being held fixed throughout the simulation. An “input” constitutes a specific Vmem pattern of the input cells. For example, input cell-1 may be set to a hyperpolarized Vmem and cell-2 to a depolarized Vmem. The input cells were then randomly set to a different state, drawn from the training batch, without disturbing the state of the rest of the network (a method referred to as “continuous computation”101) and the simulation was continued for another 300 time-steps; this was repeated for ten times. At the end of each of 10 simulations, the observed states of the output cell during the last 10 time-steps were matched against the desired output corresponding to the input in the training batch (see Figs. 2b and 3b for examples). An error was then calculated based on the difference for each of the 10 simulations for each input sequence in the batch and averaged over. The derivatives of the final average error (gradient) with respect to all the parameters were then calculated. Due to the recurrent nature of the dynamics (Fig. 1), the simulation takes many more steps than the number of layers in the network before the output reaches an equilibrium. Thus, the gradients would in principle have to computed over all those steps—a method known as “backpropagation through time” (BPTT). However, due to computational constraints and known problems like “vanishing gradient”, the gradients are typically computed over fewer time steps—this version of backpropagation is known as “truncated” BPTT.

In this work, the gradients were computed over only the last 100 steps of a simulation. The quantum of changes for all the parameters were then calculated from the gradients after scaling it with a “learning rate” that set the pace of learning. The weights and bias were assigned separate “dynamic” learning rates, where the rates were fixed in such a way that the maximum change of a weight parameter was 0.1 and the maximum change of a bias parameter was 0.01. We adopted this strategy mainly to mitigate the vanishing gradient problem, and the assumption that small changes in the parameters should cause a smooth change in behavior. We also employed a “momentum” factor, a method often used in backpropagation, that determines the extent to which the previous parameter-change vector contributed to the following vector. We used a momentum of 0.5 for the first 100 training iterations, and 0.9 for the rest; the allowed range is [0, 1]. Finally, the parameter-updates were applied (thus error was backpropagated), and the whole process was repeated for about 600,000 “epochs”. The weights were allowed to change infinitely in both directions, whereas the biases were limited within the range [0, 1]. We deemed a network as successfully trained if it achieved an average performance error of 0.0002 or below at some point during the training.

We used the Python software package called Pytorch for computing derivatives during backpropagation. This package uses a method called “automatic differentiation” to compute derivatives based on “computational graphs”—a graph that keeps track of every computation performed during a simulation, which is then swept backwards for computing the derivatives, essentially constituting an algorithm for the “chain rule” of differential calculus.

A population of 50 BEN “genomes” was instantiated. A genome consists of a combination of the network’s parameters, namely the weight matrix and the bias vector, that define a network. Henceforth, by “genome” we refer to a BEN network, for clarity. Each network consists of the following structure: an input layer consisting of 2 nodes, an intermediate “tissue” layer consisting of 25 nodes, and the output layer consisting of a single node. Initially the tissue layer was designed to be a two-dimensional lattice. Next, the tissue layer of each individual in the population was randomly rewired with a rewiring probability chosen uniformly from the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[0,1]$$\end{document}[0,1]. Thus, the population consisted of the full spectrum of tissue layer structures ranging from lattice to random. Next, the input layer and the output layers of each individual in the population were respectively connected to approximately 50% of the nodes in the tissue layer; this gave the GA a chance to explore both denser and sparser inter-layer connectivity. Every individual in the population was then assigned randomized weight and bias parameters chosen in the same range as backpropagation. They were each then backpropagated for 80 iterations in the same fashion as described above for the elementary gates, with an important difference: the inputs for these gates were transient, as they are set for the first time-step only. The final average error of each individual became its fitness score. The GA then goes through the conventional steps of selection, mutation, creating a new generation of individuals.

We used a simple flavor of GA known as the “microbial genetic algorithm”133 that capitalizes on the notion of horizontal gene transfer observed in microbes. It is essentially a form of tournament selection where two individuals in a population are picked and pitted against each other (their fitness scores are compared). We used a “geographical selection” method where individuals are placed on a 1D ring, and only geographically close individuals are picked for the tournament. The size of selection-neighborhood is known as the “deme size”, which was set to 10% in this work. The genome of the winner is transmitted to the loser at a certain “recombination rate”, which was set to 0.1. The genome of the loser is then mutated at a certain “mutation rate” (set to 0.05) and reinserted into the population. This process is repeated for 1000 generations. The individual with the best fitness score in the final generation was deemed as the best individual—the most successfully trained tissue-level logic gate.