From Physics to Pattern: Uncovering Pattern Formation in Tissue Electrophysiology

Increasing evidence points to a role for complex physical phenomena, including mechanical forces and bioelectricity, as drivers of patterning in development and regeneration. We developed a genetic algorithm-based approach to search the space of biophysical simulations for pattern-forming processes and use it to demonstrate that Turing-like patterns can arise purely bioelectrically, without requiring any variation in gene expression. We also identify several bioelectric components that can reinforce and enhance such patterns manifested in cell transmembrane voltages.


Introduction
From Turing through present day of gene-centric developmental biology, there has been tremendous success in elucidating the "chemical basis of morphogenesis" (Turing 1952)processes facilitated by gene regulatory pathways, maternal chemical pre-patterns, morphogen diffusion, and cell-to-cell interactions. Despite its earlier beginnings (Thompson 1945, Gurwitsch 1944, Burr and Northrop 1935, progress has been much slower on the physical basis of morphogenesis: Not merely regulatory phenomena combined with communication through contact and diffusion, but complex physics including electrophysiology and mechanical forces (Mammoto andIngber 2010, Beloussov 2008). It is becoming increasingly clear that development is not merely a process of unfolding, pre-programmed, chemically-encoded computation that directs spatial orchestration. Evolution has exploited a range of rich physical phenomena to coordinate individual cell activity towards the creation and repair of complex, physically active, large-scale anatomy.
One of the most exciting recent additions to the toolkit is developmental bioelectricity: communication and control of cell behavior through ion fluxes and voltage gradients (Funk 2013). It is now clear that spatio-temporal patterns of membrane resting potential across tissues regulate subsequent structure and function (Levin and Martyniuk 2017, Sullivan, Emmons-Bell, and Levin 2016, Pai et al. 2012, as well as underlie a number of patterning disorders such as birth defect syndromes (Masotti et al. 2015, Kortum et al. 2015 and cancer (Klumpp et al. 2016, Litan and Langhans 2015, Bates 2015, Barghouth, Thiruvalluvan, and Oviedo 2015. It is thus crucial to understand the origin of these patterns, the capabilities of such circuits, and ways in which they can be efficiently manipulated towards therapeutically desired goal states. Recent efforts have begun to model the fascinating relationships between bioelectric and biochemical signaling (Pietak and Levin 2016, Cervera, Meseguer, and Mafe 2016, Cervera, Manzanares, and Mafe 2015, McNamara et al. 2016. However, it's not always clear a priori which physical effects may be important and what they may be capable of. In the case of tissue electrophysiology, we have well-defined, highly conserved, modular componentschannels, pumps, and cellto-cell gap junctionswith well-known, circuit-like effects on cell membrane voltage (albeit nonlinear and time-dependent) but other interactions can potentially matter too, such as gap junction gating, chemical modulation, electrophoretic transport, changing ion concentrations, perhaps extracellular electric fields. These may lead to subtle and unanticipated phenomena. How can this be investigated?
We would like to understand what these ubiquitous, preexisting components can do on their own, without requiring the assistance of gene regulation or additional chemical morphogensmechanisms already well-known to be capable of evolving patterning phenomena. Analytical models can be of some assistance, but they quickly grow unwieldy as more physics is incorporated. To be successful, they require prior insight into how the system may be simplified for study. Alternatively, rich computational simulations incorporating a great deal of physics are available, but with dozens to hundreds of cell and tissue-dependent parameters, it's still a challenge to know where to look in order to demonstrate what patterning is possibleand harder still to convincingly demonstrate what is not possible.
In this work we demonstrate a new approach to investigating the patterning capabilities offered by a biophysical phenomenon. By combining a rich simulator with a simple genetic algorithm and fitness functions tuned to seek out elementary patterning phenomena, we're able to screen a library of bioelectric components for pattern-forming abilities. This has led to the discovery that Turing-like patterns can potentially arise through a purely bioelectric process in nonneural tissue and has identified a number of additional, voltage-sensitive components that can assist in such patterning, as well as expanding on prior results in non-neural bioelectric memory. Ultimately, we hope to solve the difficult inverse problem of automatically mapping from desired phenotype to specific circuit components and parameters (Lobo and Levin 2015) an important challenge in biomedicine and synthetic bioengineering.

Background
All healthy cells exhibit a difference in electric potential over the cell membrane ("membrane voltage", or Vmem), as a result of the action of ATP-driven transmembrane ion pumps that maintain large ion concentration gradients, along with a variety of specialized return channels through which the ions leak back through, driven by the competition between concentration gradients and voltage differences. The conductivity of these channels can be internally modulated by membrane voltage itself, as well as by chemical ligands and other signals. The resulting complex bioelectric phenomena are most well known in "excitable" neural and muscle cells, where transient opening of voltage-sensitive sodium channels leads to a rapid depolarization that propagates along the membrane. But, the basic mechanisms are present in all cells, even if lacking the characteristic combinations of channels that lead to excitability and consequent complex behavior in transient voltage excursions. Indeed, the exploitation of bioelectric circuits for information processing was discovered by evolution long before its speed-optimization as nervous systems, and even before multicellularity (Humphries et al. 2017, Prindle et al. 2015, Kralj et al. 2011. The steady-state membrane voltage is affected by a variety of cellular processes, and in turn affects other cellular processes, as well as the voltage of neighboring cells. This resting potential is known to modulate calcium influx (Deisseroth et al. 2004), MAPK signaling (Zhou et al. 2015), and the transport of charged substances, among other effects.
Cells are selectively coupled to their neighbors via gap junctions, which provide a direct cytoplasmic bridge through which ions and small molecules can flow (Mathews and Levin 2017). These gap junctions, when open, allow cells to exercise a direct influence on the resting potentials of their neighbors. A wide variety of gap junctions proteins exist, selective in permeability for different small molecules, and regulated variously by trans-junction voltage, chemical ligands, and other signals (Harris and Locke 2009).

Investigative approach
Given the emerging significance of spatial pre-patterns of membrane voltagepatterns that anticipate subsequent morphogenesiswe naturally wish to understand how they might arise. Are they simply downstream effects of traditional chemistry-and regulation-centric patterning? Or are there processes by which the natural physics of tissue bioelectricity can create these patterns entirely on their own? Some initial insight came through the development of the BETSE simulator (Pietak and Levin 2016), a sophisticated, publicly available numerical simulator specialized for tissue (non-neural) bioelectricity. BETSE models connected clusters of cells, incorporating a wide variety of detail not found in common bioelectric circuit models: separate treatment of each ion species, nonlinear GHK transport, a global electrical model, and, notably, voltage-gated modulation of gap junctions and electrophoretic transport of biologically active small molecules, among other features. BETSE can also model experimental interventions, such as surgical cuts, channel blocks, and alterations to membrane permeability. A BETSE configuration specifies channels, pumps, and gap junctions, tissue geometry, intracellular and extracellular ion concentrations, chemical ligands (if any), continuum parameters like diffusion rates and ion leakage rates (via channels and mechanisms left unspecified), and also experimental interventions. BETSE's output includes a map of transmembrane voltage across the cell cluster over time, among other variables.
Hand-crafted configurations run in BETSE demonstrated intriguing effects involving membrane voltage, including spontaneous axis induction (i.e. a gradient of polarization) and regeneration of severed gradients Levin 2016, 2017). These results were tantalizingly suggestive, but what is bioelectric patterning really capable of? What kinds of bioelectric components does it take to produce a patternwhether in evolution's hands or our own? Rather than handcrafting endless BESTE configurations, we developed a genetic algorithm to search automatically, in parallel, for configurations that would yield patterning phenomena of interest. The genetic algorithm, termed GABEE (Genetic Algorithm for Bio-Electric Exploration) (Brodsky 2017), follows a generic, "fill in the blanks" approach to evolving a template configuration file. Starting from a template, GABEE alters the combinations of channels and parameters to produce a population of distinct "individuals", which are simulated in parallel on a compute cluster and then evaluated based on the Vmem patterns that they produce. "Interesting" individuals are preferentially retained, mutated, and evaluated again.
A key challenge here, however, is defining what is meant by "interesting". Novelty selection (Lehman and Stanley 2011), aiming to map out the space of possible spontaneous patterns, appears intriguing but seems difficult to apply, given the natural run-to-run variation in simulation results and the need to avoid selecting for numerical artifacts. Alternatively, prior work searching for and classifying patterns in cellular automata often used entropy methods (Wuensche 1999, Suzudo 2004), but these are not a natural fit for noisy, realvalued physical variables, particularly so when computational limitations ensure simulations are far too small to estimate the underlying spatial probability distributions.
The simplest analogue for entropy selection might instead be selecting for the amplitude of voltage variationfor example, the standard deviation of Vmem across all cells. This indeed encourages the development of interesting spontaneous patterns, but it also turns up also a variety of artifacts: numerically unstable simulations, boundary effects, and transient behaviors such as progressive depolarization caught partway through the act. Amplitude selection also gives only a limited view of the space of possible patterns, being completely insensitive to spatial structure.
Instead, we use a set of several, more complex fitness functions, each designed to incorporate three elements:  Amplitude selection  Penalties for pathological behavior  Preference for some particular spatial (or spatiotemporal) structure These fitness functions are all hand-coded in a simple, Python-based combinator language provided by GABEE.
For scoring spontaneous patterning, we evaluate the amplitude and structure at the very end of the simulation. In this case, the most important penalty is simply an aggressive cost applied to dVmem/dt as the simulation ends (measured in RMS). This encourages stable patterns, avoiding transients and most instabilities. We also penalize voltages outside of reasonable physiological ranges, regardless of when in the simulation they appear; this catches both physiologically unreasonable configurations and additional instabilities. Together, these penalties lead to substantially cleaner results.
We can then evaluate the capabilities of a bioelectric system in terms of which types of structures it can successfully produce and how effectively. In this study, we examine three kinds of structure: spots, stripes, and bistable memory. The spot and stripe assays are both aimed at identifying spontaneous Turing-like patterns, while the memory assay tests whether the tissue can remember a patch of cells that has been experimentally depolarized and retain this state through the end of the simulation. The memory assay can be scored fairly easily, by measuring how well the cells match a uniform initial state prior to the depolarizing intervention and how well they match the intervention stencil long after the intervention is removed. The spot and stripe assays, on the other hand, are more subtle, since they must accept a wide range of pattern shapes.
For the spot and stripe assays, we would like to select for well-formed patterns indicative of generalized pattern forming capabilities, and not degenerate special cases. Pattern features no larger than a single cell are dubious and do not clearly demonstrate coordinated patterning across multiple cells. On the other hand, features as large as the entire cell cluster are ambiguous in their identity and may even reflect the effects of boundary conditions, rather than any intrinsic process. Thus, we want to select for a particular band of wavelengths, not too short, not too long. We also need to run simulations that are large enough such that spots, stripes, and different wavelengths can be distinguished.
Explicitly spectral methods such as filter banks and Fourier transforms are somewhat complicated by the hexagonal mesh, irregular boundaries, and small domain sizes found in typical BETSE simulations. Other classical image processing techniques are applicable, however. We found that total variation, defined as ∫ |∇V| dAa generalized measure of perimeter (and hence selective for short wavelengths)applied to a moving average filter over cells (selective for long wavelengths), was an effective amplitude selection measure favoring mid-range wavelengths.
The final key element is distinguishing between spots and stripes. For simple patterning systems, the skewness of the amplitude distribution is one easy measure for distinguishing spots and stripes (Shen and Jung 2005): stripes have a symmetric distribution of amplitudes, and spots have a highly skewed distribution. To select for stripes, we penalize skewness, and to select for spots, we favor absolute skewness, albeit fed through a saturation curve to discourage it from becoming too extreme.
Because of the computational heft of the BETSE simulation (a few CPU-minutes per run) and thus long GA generation times, iteratively developing and testing the fitness functions was a challenge. To ease the design process, we substituted a much simpler simulation using the Swift-Hohenberg patternforming model system (Rabinovich, Ezersky, and Weidman 2000) as a fast and easy test case until the fitness functions behaved satisfactorily. We settled on a stripes assay with a moving average radius of approximately one cell diameter and a spots assay with radius of two cell diameters (and also excluded the boundary cells).
Interestingly, it turns out that the spot and stripe assays are not equal in difficulty; given appropriate initial conditions, rough spot-like patterns can be produced from a broader class of systems than stripes, including those lacking a strong mechanism for lateral inhibition. Random snow-like patterns, possibly with some progressive coarsening or filtering, can sometimes come out looking like passable, albeit irregular spots. The stripes assay appears to be more stringent, at least as behaved on the systems tested here, with satisfactory, numerically stable solutions emerging only through local activation, lateral inhibition mechanisms.

Methodology details
The simulation template used in the experiments here includes the following ion channels with adjustable levels of expression: Passive channels -Na + and K + membrane leaks Voltage-gated channels - Nav1.6 (persistent NaV, "NaP") 1  Kv1. Ligand-gated channels - HCN2 (sensitive to cAMP) 3  rod-type CNG (sensitive to cGMP)  a hypothetical olfactory-like, cAMP-sensitive CNG 4 The following parameters are also adjustable: Chemical kineticsindependent cAMP and cGMP production/decay kinetics & gap junction diffusion rates Gap Junctions -GJ ion permeability and voltage gating sensitivity (threshold and maximum closure) BETSE is configured with a 200μm world size and a regular hexagonal cell lattice. Simulations are initialized as per the standard BETSE initialization procedure, starting with uniform membrane voltage and chemical concentrations. BETSE's static and dynamic noise features are enabled (parameter values 5.0 and 10 -7 , respectively) to inject randomness into the simulations and to reduce the influence of boundary effects. BETSE defaults are used for the remaining model parameters. BETSE is configured with a simulation time step of 0.5ms (necessary to avoid instability with Na + channels) and a well-mixed extracellular environment. The slow chemical timescales of cAMP and cGMP are artificially accelerated by several orders of magnitude in order to make simultaneous simulation with voltage-gated channels computationally feasible. Analytical modeling suggests this should not substantially impact the steady-state results (Brodsky 2018). Simulations are run for 2 seconds of simulated time, which corresponds roughly to several hours of real time taking into account the accelerated timescale.
GA runs use populations of 89 individuals with 3tournament selection. Mutations of real-valued parameters are normally distributed with a standard deviation of 10% of the total range for linearly scaled traits and 20% for logarithmically scaled traits. Crossover did not significantly improve performance and so was left disabled, leaving an asexual GA. Relatively high mutation rates (expected 3 per individual) and long evolutionary trajectories (400 generations) are chosen so that most trajectories approach a steady state and most with a template capable of solving the assays here do indeed find plausible solutions.
"Knockout" derivative templates are constructed based on the master template by disabling one or more bioelectric components. Their parameters remain as neutral degrees of freedom. Knockouts are chosen in attempt to narrow down which bioelectric components are necessary for good fitness and for particular visually identifiable properties. The different templates are separately initialized and run several times over for each assay, producing small collections of final best individuals.
The fitness functions are given as follows (detailed language documentation is available with GABEE (Brodsky 2017)):

Memory:
For each assay, abject failure is usually easy to distinguish consistently low scores and a lack of any plausible solutions under visual inspection. A somewhat arbitrary cutoff score could be picked to denote "success". However, solutions often seem to cluster in gradations of quality and not merely frequency of success, so relative comparisons based on score are of interest. By contrast, frequency of success reaching some threshold or number of generations to success seem to be rather noisy signals. Sometimes, templates with fewer adjustable components seem to have higher frequency of successperhaps because there are fewer interesting local minima to get trapped in. To make meaningful sense of the knockout experiments, we would prefer to use only measures that monotonically improve as the number of adjustable parameters increases.
Ideally, we would measure "best achievable fitness" for each template. If the GA is well-behaved, best achievable fitness should behave monotonically. However, such a maximum is difficult to determine accurately without large numbers of runs. Instead, one might compare by mean fitness, median fitness, or top quartile fitness. Given the skewed distribution of scores and small sample count, common parametric statistical tests are not applicable for comparing different templates, but a nonparametric test such as Mann-Whitney can be used to substantiate observed differences. To compare different templates, we score individual GA runs by the average of each generation's best individual fitness over the final 10 generations and then compare different populations of GA runs with Mann-Whitney. (Care must be taken to use an exact Mann-Whitney test (Marx et al. 2016) rather than an approximation valid only in the limit of large populations.) Suspicious GA runs are spot-checked by rerunning all topscoring individuals with a new random seed and a halved BETSE time step. In general, the fitness score based on the final 10 generations falls in most runs by up to 30%, reflecting natural variability and selection bias in having picked the best results. Robust, stable solutions may fall only a few percent, while numerically flimsy results fall more. GA runs where the score falls by more than a factor of 10 are discarded as spurious, likely indicating strong selection for numerical artifacts. Few runs are found to fall by intermediate amounts.

Example runs
In this section we consider a few sample runs of the GA, illustrating several kinds of bioelectric patterns it finds. Figure  1 shows three different examples evolved using the spots assay with the complete bioelectric template, while Figure 2 shows examples from all three assays, evolved under different templates. As the spots and stripes assays show, the bioelectric physics is indeed quite capable of forming intricate spontaneous patterns, exhibiting repetitive features at a characteristic scale, without requiring any involvement of gene regulatory networks. The memory assay also shows that the physics can be made to remember a pattern imposed by external means (in this case, a circle). The patterns show sharp bipartite regionalization, distinguishing cells through membrane voltage and the concentration of diffusible ligands completely independent of the cell's transcriptional state.
The leftmost example in Figure 1 shows a typical, highscoring spot pattern, while the other two examples show some less common results, which also happen to be less favored by the fitness function (convenient, though not always the case). The fitness trajectories show that a reasonable level of convergence was attained within the allotted number of generations, although there certainly are occasional exceptions.
All three trajectories show a nontrivial rate of "failure"where the numerical simulation failed to complete successfully for some individuals. The left two examples are fairly typical for failure rate (for NaV-containing trajectories), while the rightmost example is unusually fragile. Such failures are usually due to numerical instability detected within BETSE itself: despite the bounds configured in the template, the parameters provided by the GA were too aggressive and caused blow-ups when simulated. Often, the situation can be remedied by a smaller time step, but this increases simulation time proportionately. Since generation time is limited by the slowest individual, adjusting the time step per individual is not helpful. Instead, a global trade-off needs to be made, selecting a time step that can successfully simulate a wideenough range in the physical parameters, without being impractically slow. In these simulations, NaV channels, with their fast transition rates, seem to be the limiting constraints on stability. In other circumstances, a common culprit is high diffusion rates. It is often observed that trajectories like the rightmost one, with high failure rates, erratically varying fitness, and unusual voltages, are dancing on the edge of numerical feasibility and may be actively selecting for numerical artifacts; such results should be cross-checked by running winning individuals under an altered time step (successfully, in this case). Figure 2 demonstrates all three assays both under the complete template and under more limited, "knockout" templates. Unsurprisingly, there exist knockouts that eliminate the original patterning capability entirely (right column)in this case, the removal of the CNG channels for spots and stripes, and the removal of both CNG and NaP for memory. An alternative mechanism for poor quality spots remains (employing NaP), while stripes are eliminated entirely. Interestingly, however, there also exist knockout templates that do not eliminate patterning but instead merely weaken it, producing recognizable yet qualitatively altered patterns of lower fitness. In the middle column, we have soft-edged spots and filamentous stripes, lacking sharp borders and showing weaker distinction between cells inside and outside the pattern regions. The memory assay remains rather resilient when stripped down, but the spot and stripe patterning mechanisms were apparently able to exploit a variety of different voltagesensitive components in order to improve pattern intensity and regionalization. This will be examined further in the next section. Simulation failure rates (red 'x' marks, right vertical axis), largely due to numerical instability, are typical in the left and middle examples but noticeably elevated in the right example, possibly due to its heavy reliance on the voltage-gated sodium channel (present in all three but more prevalent in the third by at least an order of magnitude).

Results
Here we summarize the results of approximately 150 runs of the GA under different assays and templates. The fitness scores for the runs are illustrated in Figure 3, with the complete template on the left and various knockouts of increasing magnitude following it. The spread among solutions is broad, with a mixture of both good and mediocre solutions even among the results of the complete template. There are, however, some strongly significant differences among the knockouts, pointing to key roles played by several different types of bioelectric components.
Inspection of individual solutions shows that a great many employ the NaP channel when available, and it appears to be central to good-quality memory solutions. Every knockout lacking NaP, and only knockouts lacking NaP, show a statistically significant drop in fitness on the memory assay. The knockout lacking both NaP and CNG shows consistent, complete failure. It appears that NaP is both necessary and sufficient among the variables considered here for the bestscoring type of memory solution (while CNGs also allow a lesser-quality solution, similar to spontaneous spots). The mechanism is presumably along the lines suggested by Cervera, Mafe, et. al. (Cervera, Manzanares, andMafe 2015, Cervera, Alcaraz, and, where the positive feedback of a sodium channel's voltage gatingopening in response to the depolarization, and causing depolarization by openingyet not subsequently slamming shut like classical NaV channels doleads to "negative resistance" and cellautonomous Vmem bistability.
On the other hand, this mechanism alone is clearly insufficient for good quality spots and stripes, where CNGs appear to be critical among the components included in our template. All CNG knockouts show dramatic drops in fitness on the spot and stripe assays, with performance bordering on complete failure.
Inspection of individual solutions shows that the highscoring spot and stripe results seen here employ high concentrations of one or both CNGs, often at the upper limit. The mechanism appears to be a form what we have dubbed "autoelectrophoresis"where a charged ligand gates an ion channel in such a way that the electric potential becomes more attractive to the ligand, which then electro-diffuses inward through gap junctions from neighboring cells. This increases  the ligand's concentration locally while depleting it elsewhere a form of local activation, lateral inhibition. This effect has been observed theoretically in earlier work with BETSE (Pietak and Levin 2016), and in a companion work, inspired by early results from the GA, we demonstrate with a specialized analytical model that the mechanism is indeed capable of producing rich, Turing-like patterns (Brodsky 2018). The mechanism is shown to be closely analogous to the aggregation of bacteria via chemotaxis, as captured in the classic Keller-Segel model (Hillen and Painter 2008). Here, our results indicate that a set of real, well-characterized channelsthe CNGsare physically capable of performing this patterning feat, in conjunction with their ligands cAMP and cGMP.
It is likely that other ligand-gated channels are also capable of autoelectrophoretic patterningbut not all of them. HCN2 (also gated by cAMP), which we have included in our models, does not seem to share this ability. In nearly all examples here, it is strongly down-selected. This is despite experimental evidence that it is important in development (unpublished); presumably its significance lies in considerations not captured by our model or fitness functions. Because of its low expression and the lack of plausible solutions in the absence of CNGs, we did not attempt to manipulate HCN2 in the knockout experiments, anticipating little effect.
Interestingly, the spot and stripe assays show a much more subtle trend across the other knockouts. No single component other than CNG causes a statistically significant loss of fitness when removed. However, the simultaneous removal of NaP, KIR, Kv, and voltage-gated gap junctions (VGGJs) causes a very significant loss of fitness, giving rise to "weaker" results of the sort depicted in the middle column of Figure 2. These components are quite heterogeneous -NaP is a depolarizing channel, while KIR and Kv are hyperpolarizing, and VGGJs are neither, instead affecting the communication between neighboring cells. An obvious commonality is that they are sensitive to voltagebut so is HCN2. These components each seem to have different ways of enhancing regionalization, and apparently any one alone is often sufficient. The components and their effects share some similarities with previously observed "contrast-enhancing" components (Pietak and Levin 2017) but are even broader. NaP presumably works through its negative resistance, tuned to a lesser level than in the bistability results. VGGJs, which in their own way have a negative resistance-like property, presumably work by providing sharp isolation between neighboring cells once the voltage difference reaches a certain threshold. It is not yet clear how the K+ channels help, or even which ones are responsible. Analytical modeling suggests they might work by providing a high effective impedance for the CNG channels (Brodsky 2018)in circuit terms, an active pull-up. These voltage-gated adjuncts may also be useful even in traditional gene-regulatory patterning mechanisms for Vmem; further investigation is needed.

Conclusion
We have demonstrated the use of a simple genetic algorithm to screen a bioelectric physical simulation for pattern forming phenomena. This also serves as a proof-of-principle for the important and difficult inverse problem of translating desired outcomes into bioelectric mechanisms. We uncovered several different kinds of patterns, including the first demonstration of Turing-like bioelectric patterns, and identified the CNG and NaP channels as examples of key drivers among our ensemble Shaded regions indicate statistically significant distinguishability from full template (leftmost column), all cases p < 0.01 one-tailed Mann-Whitney. Dashed lines indicate distinguishability from adjacent column, p < 0.05 two-tailed Mann-Whitney. of test components. We also showed that a wide variety of voltage-gated components could assist in such patterning, leading to stronger, sharply regionalized patterns.