Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine

In this set of simulations, a circular cluster with limited diffusion at the outer environmental cluster boundary was simulated to explore the effect of tight junctions (Figures 10 and 11). Initial values for concentrations and voltages in the simulations were those of the final simulation for profile B cells, Table 3. Simulations investigated (1) the general effect of TJ presence, which decreased extracellular boundary diffusion constants by 1.0 × 10−5 from free diffusion values (Figures 10B and 11); (2) the effect of no TJ by leaving extracellular boundary diffusion constants at those of free diffusion values (Figure 10A); and (3) the ability for the TEP to alter cluster characteristics by affecting the permeability of voltage-sensitive gap junctions (Figure 11C), inducing electroosmotic flows (Figure 11E), and inducing self-electrophoresis/electroosmosis of membrane-bound ion pumps and channels (Figure 11D) An additional simulation, whereby the TJ barrier was broken by removal of cells during the course of a simulation, demonstrated the role of tight junctions in creating a characteristic bioelectric signal with wounding (Figure 10C).

The mere presence of TJ, even in the absence of inhomogeneous membrane channel distribution, current, or electroosmotic flow, was seen to alter cell Vmem depending on a cell’s location with respect to the cluster boundary (Figures 10 and 11). Simulations showed a trans-boundary voltage difference of approximately 24 mV (analogous to the TEP observed in organs and organisms) spontaneously appeared across the inner and outer membranes of cluster boundary clusters when ion transport between extracellular spaces at the cluster boundary was limited by simulated TJ (Figure 11B shows the TEP as a close-up to the outer cell membranes). This trans-boundary voltage gradient did not appear when ion transport in extracellular spaces was similar to free-diffusion values for ions (Figure 10A), confirming the well-known role of TJ in establishing the TEP. As detailed in Figure 11A, tight junctions generate the TEP by maintaining a transport barrier at the cluster boundary, which acts analogously to the cell membrane to internalize charge emitted by cells, leading to the generation of a typically positive voltage in the extracellular spaces internal to the cluster at the exterior membranes. As the environmental space of apical cell membranes has zero voltage, there is a voltage difference between the apical and basal membranes of outer cells, which defines the TEP.

The wounding event transitioned the cluster into a new steady state with characteristic bioelectrical pattern. Cells local to the wounded area develop depolarized Vmem in addition to a characteristic pattern of current flow featuring current directed out of the wound, which returns back to the cluster at each side of the wound (Figure 10C). This well-known pattern of endogenous current flow and associated electric field are implicated in cellular signaling events for wound healing and the initiation of limb development (Borgens, 1984; Nuccitelli, 2003b; Zhao, 2009).

Finally, when TJ are present, the electric field generated by the TEP was predicted to influence voltage-sensitive gap junction permeability (Figure 11C) to spontaneously divide the cluster into two networks: an outer layer of GJ connected cells around the perimeter of the cluster, and the inner network of the cluster bulk, with the inner and outer layers isolated by low permeability gap junctions (Figure 11C). The TEP was also found capable of inducing electroosmotic flows with, on average, velocities of about 10 nm/s (Figure 11E) and to redistribute ion pump and channels in the membranes to generate polar cell fluxes across the apical and basal membranes (Figure 11D).

We conclude that the TEP, an important bioelectrical state characteristic of multicellular clusters, can arise spontaneously in cell clusters, simply by inhibiting diffusion from extracellular spaces of the cluster boundary. The TEP contributes to the generation of a characteristic cluster-wide bioelectric state upon wounding, and creates micro-environments with differential channel activity at the boundary and interior of the cluster.

A key question in this field is the origin of bioelectric pre-patterns. Turing and others (Turing, 1952; Cross and Hohenberg, 1993) showed that pattern can spontaneously self-organize from symmetry breaking in initially homogeneous media. Could bioelectric dynamics enable voltage pre-patterns to emerge spontaneously in physiologically homogeneous tissue? Simulation 8 explored spontaneous Vmem patterning (symmetry breaking) created by a positive feedback mechanism in the networked cell cluster. A wide range of feedback mechanisms exist in bioelectrical tissue systems on account of the chemical and Vmem sensitivity of ion channels and GJ, the non-linear relationship between channel/GJ activity and Vmem, and the ability for voltage gradients to alter concentration profiles of ion channel gating ligands via electrodiffusive transport. Thus, we reasoned that positive feedback loops could amplify small physiological differences (noise) to result in macroscopic distributions that could underlie the origin of bioelectric pre-patterns. An example feedback mechanism was found capable of generating a strong, dipolar axial Vmem gradient in a cluster of voltage-sensitive GJ networked cells with TJ (Figure 12). In its initial state, the cluster had a small Vmem asymmetry of <10 mV, due to a small increase in K+ leak channels for a small number of cells in the upper right side (Figure 12A). This small asymmetrical expression of K + leak channels also corresponded to small related asymmetries in environmental voltage (Figure 12B) and anionic ligand concentration (Figure 12C). No changes to voltage-sensitive GJ open/closed state were noted in the initial state (Figure 12D, where blue = open, white = closed).

However, the presence of an electrodiffusing cationic gating ligand, which opens a Na+ channel based on its extracellular concentration (a situation analogous to that of acetylcholine), in combination with TJ and GJ activity, generated a strong axial Vmem gradient (Figure 12).

The positive feedback loop is outlined in Figure 12E. TJ maintain an environment that internalizes charge, allowing a voltage to develop in extracellular spaces (Figure 12B), while cells with more depolarized Vmem exclude the cationic ligand to the extracellular spaces (Figure 12C). As voltage in the extracellular space is the inverse of that of intracellular spaces, the gating ligand also travels extracellularly from cells with more hyperpolarized to more depolarized cells (see currents in Figure 12B), leading to further build up of gating ligand around positive cells. Positive feedback results, as the positively charged gating ligand acts to open a Na+ ion channel, leading to cell Vmem becoming more depolarized. As voltage gradients develop between cells, GJ close, further reinforcing the development of two regions of highly distinguished Vmem (Figure 12B).

We conclude that it is possible for positive feedback mechanisms to generate strong Vmem gradients a cell collective, and that this may be one mechanism through which bioelectric pre-patterns may be generated and manipulated.

Overall, BETSE enables highly detailed and accurate modeling of complex bioelectrical signals and states. A basic validation, using experimentally derived parameters, and comparing results with experimentally obtained observations for the same system (Xenopus oocytes), showed high correspondence (<10% discrepancy) between BETSE-calculated and experimental observables (Table 2). In addition, using the Nernst equation [equation (21)] to calculate Vmem on the basis of simulated electrodiffusing “reporter dye” intra- and extracellular concentrations, showed remarkable correspondence to BETSE direct-calculated Vmem. Likewise, comparison between BETSE direct-calculated and Goldman-derived Vmem also showed excellent agreement. Multicellular simulations with TJ (Figure 10) also showed the characteristic TEP across the exterior cluster boundary, with typically observed polarity (inside positive) and magnitudes (~25 mV) close to those observed experimentally (Hay and Geddes, 1985). BETSE also correctly predicted the characteristic bioelectric signal occurring with wounding (Zhao, 2009) (Figure 10D). Simulated endogenous current flows with magnitudes from 1 to 200 μA/cm2 were within the range of typically observed experimental currents 1–500 μA/cm2, Nuccitelli (1992).

Dynamical system theory maintains the concept of an attractor as a parameter (or set of parameters) toward which a system tends to spontaneously evolve, even under a wide range of starting conditions (Milnor, 1985). An attractor state is also characterized by the system’s return to the attractor state after perturbation. Simulations 2 and 3 from our validation study emphasize how a cell’s resting Vmem is an attractor state that can be reached even with remarkably variable initial conditions, such as intracellular ion concentrations being equal to those of the extracellular environment, and starting voltages being zero (Figure 5; Table 3). Also consistent with the concept of an attractor state, even with transient perturbations, such as the introduction of a high K+ concentration to the environment, cell Vmem returns to its original resting value once the perturbation ceases (Figure 6). These properties are a powerful feature of bioelectric circuits, which implement robust and stable control elements during pattern formation under a wide range of conditions (McCaig et al., 2005; Levin, 2012).

Our validation simulations also demonstrate the expected timing and shape for excitable signals, which may become possible when cells express voltage-gated Na+ and K+ channels (Figure 7). Simulations with excitable channels also highlight the important influence resting Vmem has on cell excitability, with more depolarized resting potentials being able to induce periodic self-excitations (Figure 7), which is consistent with the mechanism of pacemaker cells of the heart. Some “somatic” cells, such as embryonic amphibian epithelium and some cancer cells, are known to express voltage-gated Na+ and K+ channels; due to their lower resting potential, action-potential-like signals, including self-excitations, are predicted by our model, and have also been observed experimentally (Roberts and Stirling, 1971; Onganer et al., 2005). Moreover, some aggressive metastatic cancers exhibit a depolarized resting Vmem with abnormal, high expression of voltage-gated Na+ and a deficient expression of voltage-gated K+ channels (Onganer et al., 2005; Djamgoz, 2014). Our results (Figure 7C) also indicate that low expression of voltage-gated K+, in combination with expression of voltage-gated Na+, may permanently alter resting Vmem from a hyperpolarized (−57 mV) to depolarized resting Vmem state (−14 mV) after a single transient depolarization event, which is consistent with the abnormal resting Vmem observed in some cancers (−20 to −5 mV) (Binggeli and Weinstein, 1986; Djamgoz, 2014). A similar form of resting Vmem bistability, also developing with the action of voltage-gated channels, has been described in the bioelectrical models of Cervera et al. (2014) and Law and Levin (2015).

Our studies also highlight the significant physiological impacts that resting Vmem can exert. One-way Vmem can exert an influence on downstream patterning mechanisms is by altering the spatial distribution of important chemical signaling molecules (such as Ca2+, serotonin, or glutamate) in relation to cell Vmem, even when the compound is present at a homogeneous concentration in the extracellular bathing medium. This non-intuitive process happens because in an electrolyte medium where voltage gradients are present, ions are influenced, not only by concentration gradients (regular diffusion), but also by voltage gradients; under the right conditions they can passively move up concentration gradients. Thus, differences in Vmem alone can influence the cytosolic concentration of important signaling molecules, leading to differences in critical process, such as gene expression and enzyme function (Levin et al., 2006; Tseng et al., 2011). Our results indicate that the presence of regional Vmem differences can passively induce differences in cytosolic Ca2+ levels (Figure 8B) in a manner inversely related to cell Vmem. As Ca2+ is an important secondary messenger, and intracellular Ca2+ levels are involved in calcium-induced-calcium-release (CICR) and ion channel gating (e.g., Ca2+ gated K+ channels), these differences in the cytosolic Ca2+ concentration with Vmem may produce significant downstream effects on cell state (including enzyme function, gene expression, and apoptosis) and further evolution of bioelectric pattern (via effects on Ca2+ gated ion channels).

Similarly, effects on the concentration profile of a negatively charged signaling molecule were in direct correspondence to Vmem values (inverse to that of cationic Ca2+), which further illustrates how cell state can be influenced in divergent ways as by simply exhibiting different values of resting Vmem signaling molecules with different charge have opposite changes to their concentration profiles. Further Vmem related physiological changes, such as osmotic pressure gradients (Figure 8), can lead to differential changes in cell volume and the development of physical forces in the collective.

Due to the well demonstrated importance of cell resting Vmem states (McCaig et al., 2005; Tseng and Levin, 2013; Levin, 2014), clear comprehension of the factors involved in specifying cell resting Vmem is an essential first-step for understanding the mechanisms underlying bioelectrically mediated pattern formation and regulation. In single or isolated cells populating clusters without GJ or TJ, resting Vmem is primarily determined by the plasma membrane’s permeability to specific ions, to levels of extracellular ion concentration, and to the presence and activity of different ion pumps, such as H/K-ATPase (Veech et al., 1995; Lodish et al., 2000; Wright, 2004). BETSE accurately reproduces expected Vmem changes corresponding to these well-known factors of influence (Figures 5 and 6), with isolated cells showing strong depolarizations with increased Na+ membrane permeability, hyperpolarization with increased K+ permeability, depolarization with Cl− membrane permeability, depolarization with increased extracellular K+, and hyperpolarization with H/K-ATPase pumps and moderate to high K+ membrane permeabilities (data not shown).

However, our studies also indicate that in cell clusters where bioelectric circuits are created via the presence of GJ (enabling intercellular communication) and TJ (blocking extracellular transport at the external boundary), factors influencing the resting Vmem attractor state can be quite different from those for isolated cells. Gap junction conductivity, TJ restriction at the boundary, and the overall geometry of the cluster are additional influences of the resting Vmem state that are unique to multicellular clusters. Simulations indicate that the degree of GJ connectivity has a strong affect on the ability for cells with altered membrane permeability properties (e.g., cells with a different ion channel expression profile or open ion channel state) to manifest different resting Vmem states in the collective (Figure 9). Since research indicates that resting Vmem is a key instructive signal (Levin and Stephenson, 2012; Pai et al., 2012), GJ are indicated as important elements determining how potent differential membrane ion channel states will be in affecting physiological outcomes. Our results are consistent with recent reports of Cervera et al. (2016) and Cervera et al. (2015), whose equivalent electrical circuit model of GJ connected cells demonstrates the wide range of cell Vmem that can result from differing GJ connectivity in a collective. It is also apparent that TJ (in conjunction with GJ) are responsible for the spontaneous emergence of a TEP, whereby cells at the outer boundary of a cell cluster acquire a depolarization compared with interior cells, thereby naturally forming a spatially dependent Vmem pattern.

The involvement of GJ- and TJ-related factors in the specification of multicellular Vmem attractor states sheds light on the importance of GJ and TJ in embryonic development, and help explain why loss of GJ permeability (Leithe et al., 2006) and increase in TJ permeability (Soler et al., 1999) are associated with cancer progression. With regard to embryonic development, our results suggest functional GJ and TJ play important roles in establishing primary patterns of Vmem and asymmetry in developing clusters (Figures 9 and 10). In cancer, which is characterized by cells with depolarized Vmem, loss of GJ communication would allow cells with different membrane channel states to express dramatically altered Vmem compared to cells with high GJ connectivity (Figure 9), a conclusion that is also clearly shown by Cervera et al. (2016). Likewise, an increase in TJ permeability is predicted to result in decrease of the TEP with consequential depolarization of interior cells and atypical channel functions for cells on the boundary and interior (Figure 10). For the first time, BETSE allows quantitative study of tissue-level electric fields interacting with resting potential gradients – two key areas of developmental bioelectricity that have heretofore been studied separately.

Bioelectric circuits are rich with opportunities for feedback cycles, implying complex dynamics can exist in networks of connected cells, forming a layer of control with its own intrinsic behavior and self-organizing capabilities. Feedbacks exist because of non-linearities in the bioelectric system on a hierarchy of scales. On a fundamental level, voltages are created by net imbalances in ionic charge density; however, ionic concentrations are influenced by voltage gradients, thereby creating a primary non-linearity in the electrolytic system. Moving to the cellular and multicellular scale, the function of ion pumps and channels alters Vmem, which in turn can affect voltage-sensitive channels and electrical synapses. Furthermore, the presence of ionic messengers, including Ca2+ or anionic serotonin or glutamate, are charged molecules subject to movement in electrochemical gradients, but can also alter Vmem directly via their effect on specific ion channels as gating ligands (Levin et al., 2006; Berridge, 2014). An example of a feedback mechanism capable of inducing a strong, axial Vmem gradient in a cell cluster was demonstrated by implementing dynamics similar to those of acetylcholine (Figure 12). It will be crucial to perform quantitative simulations of specific systems to understand the origin of the instructive bioelectric pre-patterns that regulate, for example, the formation of the vertebrate face (Vandenberg et al., 2011), and anterior–posterior polarity in planaria (Beane et al., 2013). As channelopathies are increasingly observed to be an important cause of birth defects (Bendahhou et al., 2003; Barel et al., 2008; Adams et al., 2016), such models will provide not only mechanistic explanations of the origin and progression of bioelectric pre-patterns but could also be used to test prospective interventions in silico for biomedical applications.

Beyond embryogenesis, such quantitative models will reveal the conditions under which aberrant physiological states become normalized or established as tumors, and help formulate strategies for suppressing and perhaps reprograming established oncogenic states (Arcangeli et al., 2009; Chernet et al., 2015).

This first version of the BETSE platform has several limitations. First, BETSE currently considers a fixed number of cell grid points and does not compute cell division. Future work will extend the BETSE model to include the ability to model cell proliferation and apoptosis, as can occur downstream of bioelectric signaling in development, regeneration, and cancer. Similarly, cell mobility, including galvanotaxic movement of cells and galvanotropism in response to endogenous, global electric fields, is another key component of bioelectric pattern regulation. BETSE currently considers a fixed lattice of cells, which cannot move with respect to the environment. Cell and cluster shape changes in response to endogenous signals, such as the global current and field density are planned for future work. Besides the plasma membrane, intracellular organelles, such as the mitochondria, endoplasmic reticulum, and nucleus (Mazzanti et al., 2001), have their own bioelectic control mechanisms (subcellular membrane ion pumps and channels, as well as transmembrane voltages) which may interface with cell- and tissue-level bioelectric mechanisms. Models of these intracellular components are currently being developed for BETSE, in order to assist in understanding the hierarchy of interacting systems and sub-systems and their role in developmental/regenerative pattern and cancer development and regulation. Finally, we are working to add transcriptional readouts of bioelectric state change, to allow BETSE to model the control of gene expression by Vmem change and to integrate these models with existing in silico gene regulatory networks that underlie pattern regulation (Geard and Willadsen, 2009; Lobo and Levin, 2015).

Due to the complexity of bioelectric tissue systems, suitable and realistic models, such as BETSE are essential for exploring the properties and feedback mechanisms involved in bioelectrical circuits, and to elucidate the mechanisms involved in creating and regulating pattern. These tools represent an important enabling step for exploiting bioelectrical signaling in synthetic bioengineering approaches that harness self-organizing and control capabilities of voltage gradients for guided self-assembly of patterning tissues in vitro; moreover, they are a core component of forthcoming modeling tools that will identify specific manipulations of biophysical state that are predicted to achieve desired system-level outcomes (anatomical and physiological state). Future work will use BETSE to explore details of complex feedback loops involved in pattern emergence and dysregulation, with a focus on developing bioelectric interventions for rational control of morphogenesis (pattern emergence), regeneration (pattern after perturbation), and cancer development and suppression (pattern dysregulation).

ML and AP worked closely together on the functional spec and overall system strategy, based on project conceived by ML. AP designed and implemented the system, including all mathematical models, coding, and data representation. AP generated data from simulations, and analyzed it together with ML. AP and ML wrote the manuscript together.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors extend their sincere thanks to Cecil Curry for his assistance with development of BETSE software. They also thank Richard Nuccitelli, College of Biological Sciences, UC Davis, for helpful discussions relating to the nature of endogenous currents and fields. This work was supported by an Allen Discovery Center award from The Paul G. Allen Frontiers Group.

The authors gratefully acknowledge support from the National Institutes of Health (AR055993, AR061988, HD81401, HD081326), the Paul G. Allen Family Foundation, the G. Harold and Leila Y. Mathers Charitable Foundation, National Science Foundation award # CBET-0939511, the W. M. KECK Foundation, and the Templeton World Charity Foundation (TWCF0089/AB55).

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fbioe.2016.00055