Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine

1. Introduction

1.1. Bioelectricity: Why Model Electrical Activity in Non-Neural Cells?

Explaining and learning to control large-scale pattern is a central unsolved problem, with implications for mitigation of birth defects, and the advancement of regenerative medicine and synthetic bioengineering. The dynamics of signals orchestrating large-scale order in vivo are a key area of research, as understanding these signals is an essential first step in developing interventions that alter anatomical outcomes. The dynamics of chemical signals and their gradients are becoming increasingly well-understood (Reingruber and Holcman, 2014; Slack, 2014; Werner et al., 2015). However, endogenous bioelectric signals represent a parallel regulatory system that exerts instructive control over large-scale growth and form. Recent work has demonstrated that ionic and bioelectrical signaling of various cell types underpins a powerful system of biological pattern control [reviewed in Nuccitelli (2003a), McCaig et al. (2005), Levin (2012, 2014), Levin and Stephenson (2012), and Tseng and Levin (2013)]. Importantly, endogenous bioelectric gradients across tissues can be a very early pre-pattern for subsequent transcriptional and morphogenetic events. For example, during craniofacial development of frogs, specific transmembrane voltage (V_mem) patterns determine the downstream shape changes and gene expression domains of the developing face (Vandenberg et al., 2011; Adams et al., 2016) and brain (Pai et al., 2015). Furthermore, experimental modulation of cell V_mem states can radically alter large-scale anatomy, for example, inducing eye formation in ectopic body areas, such as the gut, where the master eye regulator Pax6 cannot induce eyes (Pai et al., 2012), reprograming the regeneration blastemas of planaria to produce heads instead of tails (Beane et al., 2011), or rescuing normal brain patterning despite the presence of mutated neurogenesis genes, such as Notch (Pai et al., 2015).

1.2. Local and Long-Range Order in Bioelectrical Networks

On the scale of single cells, the V_mem spanning every living cell’s plasma membrane is a demonstrated regulator of key processes, such as cell proliferation (Blackiston et al., 2009), programed cell death (Boutillier et al., 1999; Wang et al., 1999), and differentiation (Ng et al., 2010), and is known to be a factor in the activation of immune cells (Bronstein-Sitton, 2004). For example, despite the action of growth factors, stem cells have been inhibited from differentiation by preventing the cells from developing a hyperpolarized V_mem (Sundelacruz et al., 2008). The bioelectric properties of single cells are fairly well-understood (Lodish et al., 2000; Wright, 2004). However, bioelectric states often regulate large-scale anatomical properties, such as axial polarity (Marsh and Beams, 1952; Beane et al., 2011), organ size (Perathoner et al., 2014) and shape (Beane et al., 2013), and induction of formation of whole appendages (Adams et al., 2007; Tseng et al., 2010). Moreover, pattern control involves long-range coordination of bioelectric states. In metastatic conversion (Morokuma et al., 2008; Blackiston et al., 2011; Lobikin et al., 2012), tumor suppression (Chernet and Levin, 2014; Chernet et al., 2015), brain size regulation (Pai et al., 2015), and head–tail polarity in planarian regeneration (Beane et al., 2011), the patterning outcome in one region of the animal is a function of the bioelectric states of both local and remote cells. Thus, it is imperative to understand not only how ion channel and pump activity controls single-cell electrical properties but also how electrical gradients self-organize, propagate, and evolve in multicellular networks. Moreover, understanding the origin of developmental order also requires that we understand how tissue-level gradients of bioelectric properties arise.

In a multicellular collective, endogenous patterns of V_mem and electric fields provide positional information and achieve long-range coordination of cell activity. As in the central nervous system, this occurs because cells in a tissue are not isolated, but are electrochemically connected (and, therefore, communicating) in several ways, including intracellular channels known as gap junctions [GJ (Goodenough and Paul, 2009)], and by ephaptic coupling created by local field potentials, which enable one cell’s V_mem activity to influence that of its neighbor’s (Zhou et al., 2012). These connections between cells create bioelectrical circuits involving long-range signal patterns through whole structures, which have been determined crucial for developing embryos (Jaffe, 1981; Hotary and Robinson, 1990; Hotary and Robertson, 1994; Shi and Borgens, 1995), normal limb development of animals (Altizer et al., 2001), healing of wounds (Nuccitelli, 1992, 2003a; McCaig et al., 2005; Zhao, 2009), and even in continuous tumor suppression in adult animals (Chernet and Levin, 2013, 2014). The ability for cells to couple and communicate makes local changes to cell V_mem relevant in terms of long-range signals capable of affecting the whole. Likewise, the inability for cells to form communication networks, for instance, due to improper expression or function of GJ connections, is observed in disease processes, such as cancer (Leithe et al., 2006; Trosko, 2007). Even briefly altering the bioelectric connectivity of a cellular network enables rewriting of an organism’s target morphology. For example, genomically normal fragments of planarian flatworms can be induced to regenerate heads with shapes and internal anatomy belonging to other extant species (Emmons-Bell et al., 2015), or changed to a two-headed form that regenerates with two heads in perpetuity, illustrating the ability to stably re-wire bioelectric circuits with permanent changes to the overall anatomy (Oviedo et al., 2010).

Another important bioelectrical signal relevant to multicellular clusters is a voltage gradient known as the trans-epithelial potential (TEP), which forms at the outer boundary of an organ or organism. The TEP is also implicated in normal developmental processes (Shi and Borgens, 1995), wound healing (Zhao, 2009), and disease processes, such as cystic fibrosis (Hay and Geddes, 1985), fungal infection (Gow and Morris, 1995), inflammation, and cancer (Soler et al., 1999). The TEP is created when multicellular structures develop impermeable tight junctions (TJ) between cells at the exterior boundary (Hay and Geddes, 1985); disruptions to this process induce electric fields that serve as guidance cues for many migratory cell types during injury response (McCaig, 1990; Zhao, 2009; Yamashita, 2013) and limb development (Borgens, 1984; Borgens et al., 1987). Understanding plasma membrane voltage gradients and transepithelial potentials, and their spatio-temporal transitions in vivo, is a key enabling step for the field of developmental bioelectricity and its applications.

1.3. Modeling: The Need for In Silico Simulation

Understanding and learning to control patterning signals requires a quantitative appreciation of their intrinsic dynamics and the way they evolve through time. Since the pioneering work of Turing (Turing, 1952; Raspopovic et al., 2014; Watanabe and Kondo, 2015), much effort has gone into mathematical modeling of the dynamics of biochemical signals and their gradients. While there are many platforms for modeling spiking activity in the brain (Bower and Beeman, 2007), there are few available frameworks for formulating predictive models of bioelectric signaling during slower processes involved in somatic cell pattern regulation (Cervera et al., 2016), and even fewer working from the more biorealistic perspective of ionic concentrations and movements, rather than an equivalent electric circuit model. Such biorealistic models are crucial if we are to develop effective interventions that target powerful bioelectric control processes. Furthermore, ion channels and GJs are themselves voltage-sensitive (Nau, 2008; Palacios-Prado and Bukauskas, 2009). This means that cell groups can implement highly non-linear behaviors and feedback loops that are too complex to predict or control by direct inspection. While recent efforts have begun to model some of the interesting behavior of these GJ-coupled dynamical systems (Cervera et al., 2014, 2015; Law and Levin, 2015), there is a need for a flexible, powerful platform to facilitate in silico experimentation and model-building, and for connecting bioelectric dynamics with other aspects of physiology, physical forces, and genetic networks. The availability of a realistic modeling system for bioelectricity will enable (1) formulation of models of specific patterning events based on realistic physiological and channel expression data, (2) design of predicted intervention strategies for inducing desired changes in electrical state and downstream patterning outcomes, and (3) investigation of the broader capabilities of non-neural bioelectrical networks for use in synthetic biology (Doursat and Sanchez, 2014; Kamm and Bashir, 2014; Mustard and Levin, 2014) and unconventional computation architectures (Adamatzky and Jones, 2011; Adamatzky et al., 2012).

As a core component of enabling the unraveling of the bioelectrical dynamics of tissues in this exciting emerging field, we have created the Bio-Electric Tissue Simulation Engine (BETSE) to quantitatively explore bioelectrical signals in networked cell collectives. BETSE integrates a diverse range of mechanisms and physiologies to enable model building and hypothesis testing at a level congruent with experimental observables, including electrodiffusion of multiple ions under chemical and electrical gradients in various contexts; consideration of concentration, charge, voltage, and current in both intra- and extracellular networks in order to capture important signals, such as tissue-wide endogenous ion currents, TEP, and local field potentials; and dynamic control of membrane permeability and gap junction state to simulate voltage and ligand-gated channels. This work is the first in a series of studies modeling specific patterning systems, and using BETSE to infer targeted modulation strategies. Here, we discuss the design of BETSE, validate BETSE’s bioelectrical modeling performance, and provide some insights into the fundamental mechanisms involved in patterning of networked multicellular clusters.

2. Materials and Methods

2.1. Model Overview

Whether working with metals, semiconductors, or the salt-water electrolyte of biological systems, voltages (electric potential energies) are created by net electrical charge. In typical electrical systems, such as metals and semiconductors, the charge carriers are electrons or the absence of electrons (holes). In electrolytes, ions from dissolved salts can develop concentration profiles generating net charge in a region of space and, therefore, create voltages. Furthermore, mass flux of ions can generate a net current, which is associated with intracellular and tissue-wide electric fields. Therefore, ions are the fundamental units of the bioelectrical system, and their concentrations, mass fluxes, and transport mechanisms are ultimately important. BETSE can consider ions relevant to most living systems: Na⁺, K⁺, Cl⁻, Ca²⁺, ${\text{HCO}}_3^{-}$, H⁺, and charged macromolecules, such as proteins (X⁻). In addition, BETSE can consider the movement of a charged biomolecule, such as a voltage reporter dye, glutamate, serotonin or inositol triphosphate (symbolized as Yⁿ⁻ or Yⁿ⁺, where n is a variable charge number) present at low concentrations and, therefore, assumed to not affect voltage directly due to its inconsequential contribution to local charge density.

Cells create and control V_mem by selectively altering ion fluxes across their membrane. Ion pumps, such as the sodium potassium pump (Na/K-ATPase), use free-energy released from ATP hydrolysis to move ions across the insulating cell membrane, creating net ionic charge density and voltage gradients inside and outside of the cell, similar to a self-charging capacitor (Veech et al., 1995). Ion channels in the plasma membrane allow charge to move under these concentration and voltage gradients, altering charge densities and thereby changing the concentration and voltage gradients to create bioelectrical signals. At its core, BETSE keeps track of ion concentrations and ion fluxes in space and time, reducing them to net charge distributions inside and outside of the cell, using these net charges to calculate voltages inside and outside of the cellular space, calculating changes to concentrations resulting from ion mass fluxes resulting from concentration/voltage gradients and by active ion pumps, and calculating endogenous currents from the net mass flux of ions. Membrane permeability to specific ions is used as a dynamic variable to simulate the action of specific ion channels (including K⁺ leak channels, calcium gated K⁺ and Cl⁻ channels, and voltage-gated Na⁺, K⁺, and Ca²⁺ channels).

The following details how electro-diffusive transport, voltage calculations, ion pumps, ion channel dynamics, voltage-sensitive GJ, and electroosmotic flows are handled in BETSE. Further details regarding BETSE’s underlying theory and implementation can be found in Supplementary Material. Table 1 summarizes key parameter and typical variable values and their units. A highly simplified schematic of the “bioelectric circuit” implemented in BETSE is shown in Figure 1.

TABLE 1

Table 1. Main model parameters and variables.

FIGURE 1

Figure 1. The fundamental “bioelectrical circuit” implemented in BETSE, shown on a simplified geometry of two triangular cells (1 and 2) surrounded by their respective extracellular spaces (3–7). Note that in BETSE, and in contrast to the simplified image shown, cells are defined from a Voronoi diagram and are polygonal with four or more membranes, and that a larger network of 10–1000 cells is considered in simulations. Each cell–extracellular junction has a capacitive component (membrane capacitance C_m), a “resistive” component (cell membrane diffusion coefficients, D_m), and a variable current source (representing the action of pumps, i_p). Transfer between two cells occurs via GJ, which are represented by a “resistive” component (D_gj). Transfer between extracellular spaces and to the environment is handled using “resistive” components (D_o). Boundary conditions at the global environmental boundary are represented by grounded voltage (V = 0) and fixed concentrations representing an open boundary with Dirichlet conditions. Self-capacitances for each cell and extracellular space are not shown.

2.2. BETSE Platform and Performance

BETSE is a finite volume method multiphysics simulation platform, uniquely specialized to work with a range of bioelectric phenomena arising in biological tissues, which are highly spatially heterogeneous by nature.

BETSE was implemented in Python 3.4, making heavy use of the scientific and engineering toolboxes Numpy, Scipy, and Matplotlib (Millman and Aivazis, 2011).

To make each time step of a simulation as quick as possible, BETSE uses matrix-based differential equation solvers, making memory one of the limitations of simulation size and extent. Simulating a square millimeter of tissue (~10,000 cells) with all features enabled (e.g., extracellular space simulation, electroosmotic fluid flow, all ion types included) uses approximately 14 Gb of RAM, and is considered the current limit of simulation size.

BETSE code is available from the public repository: http://ase.tufts.edu/biology/labs/levin/resources/software.htm

2.3. Core Mathematical Strategy

Biological tissue represents a challenging modeling scenario due to its highly heterogeneous nature, where closely spaced (~10–30 nm), membrane bound, electrolyte-filled cells are individually interacting with a small extracellular space at individual plasma membranes, and where the extracellular spaces connect with a continuous, aqueous environment at the cell cluster boundary. Individual cells are also connected internally via transmembrane channels, such as GJ, which enable passage of small molecules and ionic current between cells. To manage this involved biophysical situation, BETSE uses an irregular Voronoi diagram-based cell grid, embedded within a regular square environmental grid, to model the heterogeneous nature of tissues, while also allowing modeling of a continuous environmental space around the cell cluster (Figure 2A).

FIGURE 2

Figure 2. BETSE computations use an irregular Voronoi-based cell grid interacting with a regular square environment grid to model heterogeneous tissues (A). The cell grid is composed of cell center (“Δ”) and membrane points (“*”) (B). Membranes that lack a neighboring cell (and, therefore, interact only with the global environment) are identified using a boundary search algorithm. The environmental grid consists of regularly spaced points (“o”) which are tagged as being internal or external to the cell cluster (B). The intercellular spaces of the cell grid are assumed to be connected to the extracellular environment by fluxes or gradients (A), which use points interpolated between the cell membrane and environmental grid midpoints with a weighting function (cell membranes seen per grid square) to properly assign the mole transfer for a mass flux between cell and environment (B).

Each modeled cell in the cell grid has a center point (indicated in grid diagrams as Δ, see Figure 2), where scalar cell properties, such as concentration (c_i) and intracellular voltage (V_cell) are defined. Each cell also has a unique volume (vol_cell) and perimeter representing the cell membrane. This allows unique membrane properties, such as V_mem, to be defined for each segment of each individual cell membrane, thereby opening the possibility for study of individual cell polarizations and self-electrophoresis/electroosmosis of membrane-bound ion pumps and channels (Jaffe, 1981; McLaughlin and Poo, 1981).

Membrane-specific scalar and vector properties are defined at each membrane segment midpoint (indicated as * in Figure 2B). Each membrane segment also has normal and tangent unit vectors. The membrane midpoints of each cell interface with the central points of local environmental grid squares (red “o” in Figure 2) via a nearest-neighbor interpolation scheme. A weighting function (cell membranes seen per grid square) is used to properly assign the mole transfer for a mass flux between cell and environment, thereby conserving mass and charge of the system (see Supplementary Material for more information).

The interconnected grid systems of BETSE, which models individual cells as discrete patches, make it possible to shape the cluster into complex forms and to cut holes into the cell cluster (before or during a simulation). Holes in the tissue represent the continuous electrolyte in the region of the hole. This enables study of simple vasculature (e.g., capillaries feeding the tissue by diffusion from the environment), cysts (such as the model shown in Figure 2A), and wounding. BETSE uses bitmaps to define the shape of the cluster, cut holes, and to assign specific properties (i.e., membrane permeability) to desired regions of the modeled tissue (see Supplementary Material).

The core mathematical operators of differential equations used in BETSE are:

• gradient (∇s_j), which calculates the degree of change in the spatial property over space at grid point j

• divergence ($\nabla \cdot {\overrightarrow{F}}_j$), which measures the amount of outward flow of a vector field from each point in space, measuring the presence of a source (+ divergence) or sink (− divergence) at grid point j, and

• the Laplacian (∇²s_j = ∇ ⋅ ∇s_j), which is most intuitively expressed as the divergence of the gradient of a scalar property. When discretized, the Laplacian is a matrix, which can be inverted to give the inverse of the operation, such that if ∇²S_j = c_j then S_j = ∇⁻²c_j.

Discrete versions of gradient, divergence, and Laplacian/inverse Laplacian were defined, using standard finite difference and finite volume techniques (Schafer, 2006), on the cell and environmental grids. These core mathematical operators were then used where required in specific differential equation expressions.

The detailed features of the cell and environmental grids, the specific definition of the above mathematical operators on each of the grids, and the interaction between the cell and environmental grid models, are discussed in Supplementary Material.

2.4. Bio-Electrochemical Mass Transport

Ion transport in bioelectrical systems is influenced by gradients of both concentration (∇c_i) and voltage (∇V), with ions passively moving by a process known as electrodiffusion – a combination of regular diffusion and electrophoretic transport. In general, electrodiffusion is described by the Nernst–Planck differential equation, describing the rate of change in the concentration c_i of an ion i with charge z_i and diffusion constant D_i:

$\frac{d{c}_i}{\mathit{dt}}=\nabla \cdot \left(D_i\nabla c_i+\frac{D_i{z}_{i}q}{k_{b}T}\nabla V-\overrightarrow{u}{c}_i\right)$ (1)

here, $\overrightarrow{u}$ is a fluid flow [e.g., electroosmotic flow field (Andreev, 2013)], q is the electron charge constant, k_b is the Boltzmann constant, and T is the temperature (see Table 1).

Ions are actively transported by pumps in the cell membrane. Both passive and active transport processes generate ion fluxes (Φ_i). These combined fluxes can lead to changes in concentration and charge density, and can generate a system-wide ionic current density $\overrightarrow{J}$.

BETSE assumes passive electrodiffusive mass transport in a multicellular cluster follows three distinct pathways: (1) transmembrane, via intra- and extracellular spaces across the plasma membrane; (2) intercellular, between cellular spaces via gap junctions; and (3) extracellular, between extracellular spaces and within the global environment (Figure 3). Active transport from ion pumps is always assumed to be transmembrane. Therefore, BETSE considers the following sources of ion flux for an ion i:

• Transmembrane, from passive electrodiffusive transport resulting from gradients between the local intra- and extracellular spaces using the Goldman–Hodgkin–Katz Flux equation (GHK Flux equation), which is derived from the Nernst–Planck Differential equation for the case of electrodiffusion across a cell membrane for a non-steady-state V_mem (Bowman and Baglioni, 1984):

${\mathrm{\Phi }}_{{\mathit{mem}}_i}=\frac{z_i{V}_{\mathit{mem}}F{D}_i}{{\mathit{R\hspace{0.17em}T\hspace{0.17em}d}}_{\mathit{mem}}}\left(\frac{c_{{\mathit{cell}}_i}-c_{{\mathit{env}}_i}\mathit{exp}\left(-\frac{z_i{V}_{\mathit{mem}}F}{\mathit{RT}}\right)}{1-\mathit{exp}\left(-\frac{z_i{V}_{\mathit{mem}}\mathit{F}}{\mathit{RT}}\right)}\right)$ (2)

here, F is the Faraday constant, R is the ideal gas constant, and d_mem is the plasma membrane thickness (see Table 1). Positive fluxes are directing mass into cells.

• Transmembrane, from active transport resulting from ion pump activity. Details of how the dynamic ion pump rates (α) are calculated are given in section 2.7:

${\mathrm{\Phi }}_{{\mathit{pump}}_i}=\mathrm{\alpha }\left(c_{{\mathit{cell}}_i},c_{{\mathit{env}}_i},V_{\mathit{mem}},t\right)$ (3)

Positive fluxes are directing mass into cells.

• Intercellular, from passive electrodiffusive transport resulting from gradients between neighboring, GJ networked cells:

${\overrightarrow{\mathrm{\Phi }}}_{{\mathit{cell}}_i}=-D_i\nabla c_{{\mathit{cell}}_i}-\frac{D_i{z}_{i}q}{k_{b}T}{c}_{{\mathit{cell}}_i}\nabla V_{\mathit{cell}}+{\overrightarrow{u}}_{\mathit{cell}}{c}_{{\mathit{cell}}_i}$ (4)

• Extracellular, from passive electrodiffusive transport resulting from gradients between neighboring environmental spaces:

${\overrightarrow{\mathrm{\Phi }}}_{{\mathit{env}}_i}=-D_i\nabla c_{{\mathit{env}}_i}-\frac{D_i{z}_{i}q}{k_{b}T}{c}_{{\mathit{env}}_i}\nabla V_{\mathit{env}}+{\widehat{u}}_{\mathit{env}}{c}_{{\mathit{env}}_i}$ (5)

FIGURE 3

Figure 3. Electrodiffusive mass transport in a GJ networked cell cluster is assumed to follow three pathways (A): (1) transmembrane – between intra- and extracellular spaces across the plasma membrane; (2) inter-cellular – between cellular spaces via GJ; and (3) environmental – between extracellular spaces and in the global environment. The degree of movement of ions in both chemical and electrical gradients is handled using spatially varying diffusion coefficients, which reduce from free-diffusion coefficients in the environmental space to minimal values across simulated tight junctions and plasma membranes (A). In addition to concentration gradients, ions are assumed to move under the influence of voltage gradients [electric fields (B)]. Strong electric fields are assumed to exist on a microscopic scale across membranes and gap junctions due to heterogeneous charge distribution at membrane interfaces separated by tens of nanometers [(B), Microscopic fields]. Weaker, mesoscopic scale (10–100 μm) electric fields are assumed to be generated by net ion currents in the intra- and extracellular spaces [(B), Mesoscopic fields].

Changes in concentration are made by assuming the concentration change in an ion i depends on the divergence of the net (${\mathrm{\Phi }}_{{\mathit{tot}}_i}$) sum of all fluxes of the ion entering or changing in a particular region of the space (i.e., cells or environment):

$\frac{\partial c_i}{\mathrm{\partial t}}=-\nabla \cdot {\mathrm{\Phi }}_{{\mathit{tot}}_i}$ (6)

Net ionic charge density was calculated by summing all ion concentrations at a region of space:

${\mathrm{\rho }}_e={\displaystyle \sum _i}F{z}_i{c}_i$ (7)

The dynamics of ionic charge density were calculated from the mass flux of all ions:

$\frac{\partial {\mathrm{\rho }}_e}{\mathrm{\partial t}}={\displaystyle \sum _i}F{z}_i{\mathrm{\Phi }}_i$ (8)

The total current density of the environment or cell, $\overrightarrow{J}$, was calculated using the continuity equation in combination with the assumption of bulk electro-neutrality for electrolytes due to charge screening. Using the Continuity equation for current, the current density in a region follows:

$\nabla \cdot {\overrightarrow{J}}_o+\frac{\partial {\mathrm{\rho }}_e}{\mathrm{\partial t}}=0$ (9)

As electro-neutrality (zero net charge density) must be preserved in the bulk electrolyte, the base current density calculated by BETSE (${\overrightarrow{J}}_o$) was corrected by assuming that an internal electric field develops in the bulk electrolyte as a result of charge screening, which is the negative gradient of an electric potential φ_int:

$\overrightarrow{J}=J_o-\nabla {\mathrm{\psi }}_{\mathit{int}}$ (10)

Substituting equation (10) into equation (9) and rearranging to solve for the internal electric potential:

${\mathrm{\psi }}_{\mathit{int}}={\nabla }^{-2}\left(\nabla \cdot {\overrightarrow{J}}_o+\frac{\partial {\mathrm{\rho }}_e}{\mathrm{\partial t}}\right)$ (11)

After obtaining φ_int, it is used with equation (10) to produce the corrected current density for the system. Current density in the environment and in cell spaces was treated as separate.

Note that as movement in both concentration and electrical gradients can occur, the transport properties of bioelectrical systems cannot be strictly reduced to electrical constants, such as resistance or conductance. However, examining the Nernst–Planck equation [equation (1)] reveals that the diffusion coefficient D is able to serve as the constant of proportionality for movement in both chemical and electrical terms. In the absence of a concentration gradient, and multiplying by F z to convert mass flux to ionic current density, the Nernst–Planck Flux equation reduces to:

${\overrightarrow{J}}_i=-\frac{F{D}_i{z}^{2}q}{k_{b}T}{c}_i\nabla V$ (12)

Noting that the definition of an electric field is $\overrightarrow{E}=-\nabla V$, equation (12) parallels the equation relating current density to electric field via media conductivity $\frac{1}{\mathrm{\gamma }}$:

$\overrightarrow{J}=\frac{1}{\mathrm{\gamma }}\overrightarrow{E}$ (13)

Therefore, BETSE makes use of diffusion constants to characterize ion transport in different regions of the multicellular cluster, but can approximate conversions between conductivity and the diffusion constant.

Note that for movement across a membrane with thickness d_mem, the permeability of the membrane is simply P_mem = D_mem/d_mem.

2.5. Biological Voltages

2.5.1. Bioelectric Voltage Calculations Using a Maxwell Capacitance Matrix

The Poisson equation ($V=-\frac{{\mathrm{\rho }}_e}{\mathrm{\epsilon }}$, where ρ_e is electronic charge density and ε is medium electrical permittivity) is typically used to determine voltage from charge density. In air, a charged object will emanate a voltage gradient (electric field) into the space around it according to the Poisson equation. However, electrolytes are more complex. Due to the presence of mobile, oppositely charged ions in electrolytes, objects with steady-state voltages or bound charges collect an opposite surface charge from the electrolyte to form an electrical double layer approximately 1-nm thick in biological systems, which screens the voltage/charge of the object and a prevents long-range electrical field from developing at macroscopic distances into the electrolyte (Bazant and Squires, 2004). Moreover, biological systems are highly heterogeneous, with opposite-sign charge distributed at intra- and extracellular interfaces of the plasma membrane. This means opposite sign charges are separated by the small membrane thickness (3.5–9 nm), and that in a collective of many cells with closely interfacing membranes, charges are present in the low-volume extracellular space that is approximately 5–50-nm wide between cells. Therefore, a new technique was adopted to model voltage in the biological tissue. Voltages in the intra- and extracelluar spaces (V_intra, V_extra), and the related V_mem = V_intra − V_extra, were calculated from net ionic surface charge distributions using a formulation called the Maxwell Capacitance Matrix (Clements et al., 1975; Heinzel, 2008).

Capacitance is typically known in terms of an electrical device characterized by two metal plates (electrical conductors) with equal and opposite charge (±Q) on either plate, which are separated by a layer of insulating material (Figure 4A). The capacitance (C) is defined by the ratio of the voltage (V) between the plates in relation to the charge Q on each plate (Figure 4A):

$C=\frac{Q}{V}$ (14)

FIGURE 4

Figure 4. A capacitor is a device commonly characterized by two metal plates (conductors) with equal and opposite charge Q on either plate, separated by an insulating layer (A). A surface at a static negative voltage V, submerged in an aqueous electrolyte, does not produce a long-range electrical field (voltage gradient) in the electrolyte due to charge screening and the formation of the electrical double layer (B). The ratio of screening charge Q to surface voltage V is a self-capacitance for the system (B). Multiple insulator-separated conductors with variable amounts of charge and voltage on each conductor form a capacitive system, with the relationship between charge and voltage on the conductors described by a Maxwell Capacitance Matrix (C). Likewise, the interface between two cells can be reduced to a capacitive electrical system consisting of three conductive spaces with net charge and voltage (2 intracellular, 1 extracellular) separated by two insulators (2 cell membranes, capacitance C_mem), where each electrolyte-filled space has a self-capacitance related to electrolyte screening (D).

However, arrangements of conductors can involve capacitance via multiple insulator-separated conductors with variable amounts of charge and voltage on each conductor (Figure 4B). For the case of multiple conductors, the basic capacitance relation shown in equation (4) must be extended to parameterize more complex arrangements. For instance, for the three conductors shown in Figure 4C, the charge Q on each conductor can be related to voltages and capacitances of the system as:

$\begin{array}{l}\hfill \begin{array}{rl}{Q}_A& =C_m\left(V_A-V_o\right)+C_{\mathit{sa}}{V}_A\\ Q_B& =C_m\left(V_B-V_o\right)+C_{\mathit{sb}}{V}_B\\ Q_o& =C_m\left(V_o-V_A\right)+C_m\left(V_o-V_B\right)+C_{\mathit{so}}{V}_o\end{array}\end{array}$ (15)

here, C_m is a capacitance connecting the conductors, and C_sa, C_sb, and C_so are the self-capacitances of the conducting objects.

The self-capacitance of a conductor describes how much charge is acquired on the conductor/unit voltage applied to its surface. For electrolytes, self-capacitance is related to the ability of the electrolyte to screen voltage on a submerged object (see Figure 4B). Using Boltzmann relations for low voltages, in an electrolyte the ionic charge density ρ_e forming near the surface of an object with a surface voltage φ_o can be expressed as a function of distance, x, away from the surface as:

${\mathrm{\rho }}_e\left(x\right)=-{\mathrm{\epsilon }}_o{\mathrm{\epsilon }}_r{\mathrm{\kappa }}^2{\mathrm{\phi }}_o\mathit{exp}\left(-\mathrm{\kappa }x\right)$ (16)

where κ is the inverse Debye length for the electrolyte, which assuming the approximation for a symmetric monovalent electrolyte with total molar concentration C_tot for a typical BETSE system is expressed:

$\mathrm{\kappa }=\sqrt{\frac{2F^2{C}_{\mathit{tot}}}{{\mathrm{\epsilon }}_r{\mathrm{\epsilon }}_{o}RT}}=1.2x{10}^9\frac{1}{m}$ (17)

Integrating equation (16) with respect to x from the surface at 0 to infinity, and dividing by the surface voltage to approximate a self-capacitance/unit surface area for surfaces in the electrolyte:

$c_{\mathit{self}}=\frac{{\mathrm{\rho }}_e}{{\mathrm{\phi }}_o}=\mathrm{\epsilon }{\mathrm{\epsilon }}_r\mathrm{\kappa }=0.86\frac{F}{m^2}$ (18)

The system of linear equations derived when considering more complex arrangements of insulator-separated conductors can be expressed in a matrix form (Heinzel, 2008). For the highly simplified system shown in Figure 4B, and described by equation (15), the Maxwell Capacitance matrix (M) interrelating charge Q_k and voltage V_k on each conductor is:

$\begin{array}{rl}\left[\begin{array}{c}\hfill Q_A\hfill \\ \hfill Q_B\hfill \\ \hfill Q_o\hfill \end{array}\right]& =\left[\begin{array}{c}\hfill \begin{array}{ccc}\hfill C_m+C_{\mathit{sa}}\hfill & \hfill 0\hfill & \hfill -C_m\hfill \\ \hfill 0\hfill & \hfill C_m+C_{\mathit{sb}}\hfill & \hfill -C_m\hfill \\ \hfill -C_m\hfill & \hfill -C_m\hfill & \hfill 2C_m+C_{\mathit{so}}\hfill \end{array}\hfill \end{array}\right]\left[\begin{array}{c}\hfill V_A\hfill \\ \hfill V_B\hfill \\ \hfill V_o\hfill \end{array}\right]\\ \overline{Q}& =M\overline{V}\\ \overline{V}& =M_{\mathit{inv}}\overline{Q}\end{array}$ (19)

For the case where the set of charges $\overline{Q}$ are known, the corresponding set of voltages $\overline{V}$ can be found by calculating the pseudo-inverse of the Maxwell Capacitance matrix (M_inv) using a singular-value decomposition method.

In the biological system, we propose every point of contact between two cells represents a situation similar to the one shown in Figure 4C. In order to calculate voltage within the closely-spaced intra- and extracellular regions, and to thereby derive V_mem for a cell cluster, each cell–cell interface is reduced to a capacitive electrical system consisting of three conductors: two intracellular spaces and one extracellular space, which are each separated by two cell membranes with capacitance C_mem. Each space has net ionic charge Q_k and voltage V_k. The self-capacitance of each space is related by equation (18). This allows a Maxwell Capacitance Matrix identical to the one defined in equation (19) to be constructed for a single cell–cell junction (Figure 4D).

In BETSE, a typical cell cluster consists of many hundred to thousand cell–cell interfaces and, therefore, has a very large M, the pseudo-inverse of which was used to calculate voltage in each intra- and extracellular space from net charge Q_k in each region. To complete the calculation, V_mem are calculated by taking the difference between the intra- and extracellular voltages at a respective membrane point.

Note that the use of the Maxwell Capacitance Matrix to derive V_mem is only one component of the computation of bioelectrical variables – a simplified bioelectrical “circuit” is shown in Figure 1, and must also include electrodiffusive transport of ions via transmembrane, intercellular, and extracellular networks, in addition to active transport of ions by pumps, as described in section 2.4.

2.5.2. Assumptions Regarding the Biological Electric Field

It is important to clarify that while it is well known that the ions of electrolytes screen voltages arising from static charge distributions, thereby preventing electric field (voltage gradient) from static charge distributions from being seen past the electrical double layer, any net ionic current density ($\overrightarrow{J}$) arising from ion fluxes in the biological system is known to generate a small magnitude observable macroscopic electric field (${\overrightarrow{E}}_{\mathit{global}}$) according to:

${\overrightarrow{E}}_{\mathit{global}}=\mathrm{\gamma }\overrightarrow{J}$ (20)

where γ represents the media resistivity of approximately 0.02 Ωm. These endogenous currents and related global electric fields have been observed directly and are on the magnitude from 1 to 1000 $\frac{\mathrm{\mu }\text{A}}{{\text{cm}}^{\text{2}}}$ and exist in the extra- and intracellular spaces (De Loof, 1985; Nuccitelli, 1992, 2003a,b; Altizer et al., 2001).

BETSE assumes the existence of two types of electric field (voltage gradient) in the biological tissue (Figure 3B). At the microscopic scale (i.e., 10 s of nanometers), very strong voltage gradients are assumed to exist across membranes, gap junctions, and between extracellular spaces (especially across tight junctions) due to the presence of charge at interfaces separated by distances of 10 s of nanometers (Figure 3B). These electric fields, with strength on the order of 0.01–1.0 million volts/meter, are assumed to be the primary drivers of ion flux across membranes and junctions, however, are very short acting due to electrolyte screening. On mesoscopic scales (i.e., 10–100 s of micrometers) net ionic current density is assumed to be associated with a longer range, weaker electric field via equation (20) (Figure 3B), which is of much lower strength on the order of 0.2 volts/meter. BETSE assumes the current densities in the environment and in the cell networks are separate.

2.6. Standard Equations for Voltage Cross-Check and Validation

2.6.1. Nernst Equation

In cell physiology, two additional equations have been derived from the Nernst–Planck equation for use in specific situations involving transport across a membrane: the Nernst equation (Matthews, 2013a) and the Goldman equation (alternatively known as the GHK equation) (Matthews, 2013b).

For the case where the system consists of two compartments separated by a semi-permeable membrane, and the system is at steady-state with both zero ion flux and zero current across the membrane, the Nernst equation (21) can be used to predict the voltage or ratio of concentrations across a membrane:

$V_{\mathit{mem}}=\frac{\mathit{zR}T}{F}\mathit{ln}\left(\frac{c_{\mathit{ext}}}{c_{\mathit{int}}}\right)$ (21)

Note that the Nernst equation (21) should technically only be used for steady-state situations with no flux or current of the ion c. A suitable situation would be the equilibrium concentration of a substance such as a reporter dye, which is present at low concentrations and not subjected to active pumping by membrane transporters. BETSE uses the Nernst equation with internally computed intra- and extracellular concentrations of a passively electrodiffusing substance (i.e., modeled reporter dye) to obtain an alternative value for V_mem, which is used as a cross-check of BETSE-derived concentrations and V_mem calculations.

2.6.2. Goldman Equation

The Goldman equation applies for cases where there is a net flux of ions across the membrane, however, the net current is zero, leading to a steady-state or “resting” V_mem. Due to the action of active ion pumping in living cells, the steady-state V_mem represents a dynamic equilibrium with net ion flux but zero current, and can be estimated from the Goldman equation as:

$V_{\mathit{mem}}=\frac{RT}{F}\mathit{ln}\left(\frac{{\displaystyle \sum _i}{P}_{{\mathit{mem}}_i}{c}_{{\mathit{ext}}_i}^{+}+{\displaystyle \sum _j}{P}_{{\mathit{mem}}_j}{c}_{{\mathit{int}}_j}^{-}}{{\displaystyle \sum _i}{P}_{{\mathit{mem}}_i}{c}_{{\mathit{int}}_i}^{+}+{\displaystyle \sum _j}{P}_{{\mathit{mem}}_j}{c}_{{\mathit{ext}}_j}^{-}}\right)$ (22)

In the Goldman equation (22), ions are separated into anions (c⁻) and cations (c⁺) with concentrations inside (c_int) and outside (c_ext) of the cell membrane. The membrane has a specific permeability P_mem for each unique ion. The Goldman equation is also known as the Goldman–Hodgkin–Katz (GHK) equation.

Note that as the Goldman equation is derived from the Nernst–Planck equation (Matthews, 2013b), the Goldman equation cannot be used to accurately calculate V_mem without developing a circular dependency between concentration and voltage due to an insufficient number of degrees of freedom. Also, the Goldman equation only supplies the transmembrane voltage difference across the membrane, and does not give absolute values for the intra- and extracellular voltages, which are important for calculating cluster-wide bioelectrical signals and states, such as the TEP. However, model parameters computed in BETSE (ion concentrations and membrane permeabilities) were used with the Goldman equation (22) to cross-check and compare final V_mem values obtained using BETSE’s Maxwell Capacitance Matrix voltage solving method defined in section 2.5.

2.7. Ion Pumps

Ion pumps were modeled as enzymes using standard Michaelis–Menten enzyme kinetic relations, with reaction rates determined by thermodynamic arguments.

The equilibrium constant of a reaction, K_eqm, can be expressed both in terms of the reaction free energy under standard conditions, $△G_{\mathit{react}}^{\mathit{eqm}}$, and in terms of the reaction’s product concentrations (index k) and those of its reactants (index j) where a_k and a_j represent coefficients of stoichiometry for the reaction (Beard and Qian, 2007; Pekar, 2015):

$K_{\mathit{eqm}}=\mathit{exp}\left(-\frac{△G_{\mathit{react}}^{\mathit{eqm}}}{RT}\right)=\frac{\sum c_k^{a_k}}{\sum c_j^{a_j}}$ (23)

The electrochemical potential of a substance at concentration c_i with charge z_i in a region where there is a voltage V is expressed:

${\mathrm{\mu }}_i={\mathrm{\mu }}_o+RT\mathit{ln}\left(c_i\right)+z_{i}FV$ (24)

Furthermore, the overall free-energy of a reaction is described as the sum of the (electro)chemical potentials of its products (index k) minus those of its reactants (index j) where a_k and a_j represent coefficients of stoichiometry for the reaction:

$\mathrm{\Delta }{G}_{\mathit{reaction}}=\sum a_k{\mathrm{\mu }}_k-\sum a_j{\mathrm{\mu }}_j$ (25)

Using the Na/K-ATPase pump as an example, the overall reaction for the Na/K-ATPase pump is:

$3c_{\mathit{Na}}^{\mathit{in}}+2c_K^{\mathit{out}}+\mathit{ATP}\leftrightarrow 3c_{\mathit{Na}}^{\mathit{out}}+2c_K^{\mathit{in}}+\mathit{ADP}+P$ (26)

From the abovementioned fundamental chemical principals, the overall free energy, ΔG_pump, for the Na/K-ATPase pump reaction can be expressed (Smith and Crampin, 2004):

$\begin{array}{l}\hfill \begin{array}{rl}\mathrm{\Delta }{G}_{\mathit{pump}}& =△G_{\mathit{ATP}}^o+RT\mathit{ln}\left(\mathrm{\Omega }\right)-F{V}_{\mathit{mem}}\\ \mathrm{\Omega }& =\frac{\mathit{cADP}\mathit{cP}\mathit{cN}{a}_{\mathit{out}}^{3}c{K}_{\mathit{in}}^2}{\mathit{cATP}\mathit{cN}{a}_{\mathit{in}}^{3}c{K}_{\mathit{out}}^2}\\ \end{array}\end{array}$ (27)

when ΔG_pump = 0, the reaction is at equilibrium. Using equation (23), an expression for the Na/K-ATPase pump reaction equilibrium constant in terms of the standard free energy for ATP hydrolysis and cell V_mem is:

$K_{\mathit{NaKATP}}^{\mathit{eqm}}=\mathit{exp}\left(\frac{-\mathrm{\Delta }{G}_{\mathit{ATP}}^o+F{V}_{\mathit{mem}}}{RT}\right)$ (28)

Following with basic Michaelis–Menten enzyme kinetics, an estimate for the rate of the reversible enzymatic pump reaction follows as:

$\begin{array}{llll}\hfill {\mathrm{\alpha }}_{\mathit{NaKATP}}& ={\mathrm{\alpha }}_o\left(\frac{\frac{\mathit{cATP}}{K_{\mathit{ATP}}}\frac{{\mathit{cNa}}_{\mathit{in}}}{K_{\mathit{Na}}}\frac{c{K}_{\mathit{out}}}{K_K}}{\left(1+\frac{\mathit{cATP}}{K_{\mathit{ATP}}}\right)\left(1+\frac{{\mathit{cNa}}_{\mathit{in}}}{K_{\mathit{Na}}}\right)\left(1+\frac{c{K}_{\mathit{out}}}{K_K}\right)}\right)& \hfill & \\ \hfill & \times \left(1-\frac{\mathrm{\Omega }}{K_{\mathit{NaKATP}}^{\mathit{eqm}}}\right)\end{array}$ (29)

Values for the Michaelis constants K_Na = 5.0, K_K = 0.2, and K_ATP = 0.15 were obtained from references (Munzer et al., 1994; Vrbjar et al., 1994). Values of α_o were roughly calibrated to Na/K-ATPase pump rates reported for Xenopus oocytes (Costa et al., 1989).

In addition to Na/K-ATPase pumps, BETSE can optionally simulate Ca-ATPase, H/K-ATPase, and V-ATPase pumps using free-energy regulated pumping rates analogous to that outlined above for the Na/K-ATPase pump.

2.8. Voltage-Gated Channels

A range of voltage-gated channel types have been implemented in BETSE using Hodgkin–Huxley style differential equations to define the state of membrane diffusion to a specific ion (e.g., Na⁺) as a function of V_mem and time. Specific parameters and functional relations were obtained from the online database, Channelpedia (Ranjan et al., 2011).

The present work specifically uses a combined generic voltage-gated sodium channel (NaV) from (Hamill et al., 1991), and a delayed-rectifier voltage-gated potassium channel (KV1.2) from (Sprunger et al., 1996), to generate excitable signals. A standard Hodgkin–Huxley style model uses an electrical equivalent circuit equation to determine changes to current and voltage across a membrane, with a set of differential equations controlling the conductance of the membrane (Nelson, 2004). Since conductance is proportional to the membrane diffusion constant for a particular ion [see equations (12) and (13)], BETSE uses the same Hodgkin–Huxley style equations developed to describe membrane conductivity state to describe the membrane diffusion state of a particular ion, updating subsequent changes to currents and voltages using its own methods, as described in the above. Details regarding voltage-gated channel dynamics are specified in Supplementary Material.

2.9. Gap Junctions

Gap junctions were modeled as (optionally) voltage-sensitive conduits influencing the intercellular diffusion coefficient for all ions uniformly via a diffusion–constant scaling factor, ${\mathrm{\beta }}_{\mathit{GJ}}^o$. Simulated transport through GJ used the Nernst–Planck equation [equation (1)] to update concentration of all ions moving under intercellular concentration and voltage gradients. In the absence of GJ, cells were modeled to have an intercellular diffusion coefficient of zero (${\mathrm{\beta }}_{\mathit{GJ}}^o=0$). Medium-high GJ connectivity corresponded to ${\mathrm{\beta }}_{\mathit{GJ}}^o$ = 1.0 × 10⁻⁶, an intercellular diffusion coefficient of approximately 1.0 × 10⁻¹⁵m²/s. Assuming 1.0 × 10⁵ GJ per cell, and cylindrical GJ with pore diameter of 1.5 nm and length of 26 nm, this corresponds to individual GJ conductance of 68 pS, which is in the mid-range of reported GJ conductances (Goodenough and Paul, 2009).

Voltage gating of GJ was described using the kinetic model of (Harris et al., 1983), which calculates GJ open/closed state (β_GJ) dependence on voltage difference across the gap junction (V_GJ) and time. Specific details regarding voltage gating of GJ are described in Supplementary Material.

2.10. Tight Junctions

Multicellular organs and organisms develop very low-permeability TJ at their exterior boundary, which are involved in creating the important TEP voltage gradient across the organ/organism boundary. In BETSE, the degree of movement of ions in both chemical and electrical gradients was handled by considering three interconnected, but distinct transport pathways (transmembrane, intercellular, extracellular), with the possibility for spatially varying diffusion coefficients within extracellular regions, with low diffusion at the boundary simulating the presence of TJ (see Figure 3).

2.11. Electroosmosis

Electroosmotic flows are a hypothesized transport mechanism in biological systems (Andreev, 2013). BETSE assumed that electroosmosis may occur through small channel structures of the heterogeneous tissue, such as gap junctions between cells (gap junction radius r_gj ~ 5 to 8 nm) and the narrow channels (d_ecm ~10 to 30 nm) formed by extracellular spaces.

Our simple estimate used a modified version of the Hagen–Poiseuille equation (Gao et al., 2011) to estimate electroosmotic fluid flows between the small channels represented by gap junction connected cells or extracellular spaces:

${\overrightarrow{u}}_o=\frac{\mathrm{\pi }{r}^4}{8\mu }{\overrightarrow{F}}_e$ (30)

where ${\widehat{F}}_e$ is a volume force generated by electrostatic forces resulting from a voltage gradient (electric field $\overrightarrow{E}$) between two cells or extracellular spaces:

${\overrightarrow{F}}_e={\mathrm{\rho }}_e\overrightarrow{E}$ (31)

As mass cannot be created or destroyed, fluid flow velocity must be a divergence-free field, which physically corresponds to the development of internal pressures resisting fluid flow. The internal pressure was estimated as:

$P_{\mathit{int}}={\nabla }^{-2}\left(\nabla \cdot {\overrightarrow{u}}_o\right)$ (32)

The gradient of the internal pressure was used to correct the velocity calculated from equation (30), yielding the final estimate for electroosmotic fluid velocity:

$\overrightarrow{u}={\overrightarrow{u}}_o-\nabla P_{\mathit{int}}$ (33)

Electroosmotic fluid velocities were treated separately in the intra- and extracellular spaces.

2.12. Other Biophysical Phenomena

Details regarding the implementation of other biophysical phenomena, such as lateral self-electrophoretic/electroosmotic transport of ion pumps and channels in cell membranes, and the development of osmotic and hydrostatic pressures, are discussed in Supplementary Material.

3. Results

3.1. Model Validation and Resting V_mem Regulation in Isolated Cells

Simulations 1, 2, and 3 were used to validate the core BETSE model by determining its ability to predict resting V_mem and expected V_mem dynamics under a series of perturbations for isolated cells not connected by TJ or GJ (single-cell behavior). Validations also checked that equilibrium concentration profiles of an electrodiffusing charged molecule (simulated reporter dye) showed values predicted from the Nernst equation (21). The behavior of voltage-gated channels was explored in simulation 4.

3.1.1. Simulation 1: Prediction of Xenopus oocyte V_mem

The first validation step used experimentally derived input values (membrane permeability and environmental ion concentrations), comparing simulated output to experimentally observed parameters (V_mem).

Experimentally observed membrane ion permeabilities and extracellular ion concentrations of Na⁺, K⁺ and Cl⁻ obtained from Xenopus oocytes (Costa et al., 1989), were used as input parameters (Table 2). The simulation was performed on a small network of 35 isolated cells for 30 min of simulated time. The resulting BETSE-derived V_mem and intracellular ion concentrations were compared to those observed experimentally for Xenopus oocytes with the same membrane ion permeabilities and under the same extracellular ion concentrations (Costa et al., 1989). After 30 min of simulation, steady-state V_mem and intracellular ion concentrations calculated by BETSE showed <10% difference between experimentally determined values measured from Xenopus oocytes (Table 2).

TABLE 2

Table 2. Input membrane permeabilities (P_{mem_i}) and simulated V_mem, and ion concentrations in the extra- and intracellular spaces at steady-state (30 simulated minutes) for the resting membrane ion permeability profiles of Xenopus oocytes bathed in Ringer’s solution, as reported elsewhere (Costa et al., 1989).

3.1.2. Simulation 2: Resting V_mem as an Attractor State

Simulation 2 explored resting V_mem as a dynamic systems attractor state, reaching a characteristic value even with highly divergent initial conditions. This is an important property to understand, in light of the significant robustness of biological pattern regulation. The simulation was performed on a small cluster of 183 isolated cells, which were not connected by gap or tight junctions, where cells in different regions were assigned to one of three membrane ion permeability profiles (Figure 5A). The membranes of profile A, B, and C cells had high, medium, and low sustained K⁺ membrane permeability, respectively (Figure 5). All other parameters associated with cells in the three profiles were identical.

FIGURE 5

Figure 5. Resting potentials (steady-state V_mem) are states of dynamic equilibrium highly influenced by cell membrane ion permeability profiles. Three different cell membrane ion permeability profiles are defined for cells in a cluster (A). Starting from equal concentrations of ions in the intra- and extracellular environments (see Table 3), BETSE reaches a steady-state V_mem value closely predicted by the Goldman equation for the three cell profile types (A,B).

The simulation began with a non-physiological starting state featuring equal concentrations of ions in both the intra- and extracellular environments (starting concentrations are those typical of human plasma and are given in Table 3) and with no voltage in any part of the system (V_mem = 0 in all cells).

TABLE 3

Table 3. Simulated ion concentrations in the extracellular and intracellular spaces at time zero, and at steady-state (20 simulated minutes) for three different membrane ion permeability profiles defined in Figure 5.

The Goldman equation [equation (22)] was used with membrane permeability and ion concentration values available at each time step to predict V_mem using conventional measures and provide an indicator of expected resting V_mem for each of the three profiles.

The simulation shows that after 20 simulated minutes, the BETSE-calculated V_mem closely approaches (<10% discrepancy) the value predicted by the Goldman equation [equation (22)] for the three cell membrane-permeability profile types (Figure 5). As is expected from theory (Matthews, 2013b), the steady-state V_mem value complying with the Goldman equation is reached when the net trans-membrane current reaches zero (data not shown).

In addition to the six major ions, an electrodiffusing negatively charged “reporter dye” was also included in the simulation (“Dye⁻,” Table 3) and assumed to be at low concentrations and, therefore, not influencing V_mem. The Nernst equation [equation (21)] was used with BETSE-simulated intra- and extracellular dye concentrations as an alternative V_mem estimate (“V_mem Dye,” Table 3); results are virtually identical between the direct-BETSE and dye-estimated V_mem values.

Notably, while concentrations in intra and extracellular spaces began equal, at steady-state (20 simulated minutes) intracellular ion concentrations were within expected physiological ranges (Veech et al., 1995; Lodish et al., 2000; Wright, 2004) (Table 3).

In addition to model validation, this simulation emphasizes resting V_mem of isolated cells as stable states of dynamic equilibrium that are attractor states with final values highly influenced by cell membrane ion permeability profiles. As expected, increased membrane permeability to K⁺ (simulating increased expression of K⁺ leak channels) leads to higher degrees of V_mem hyperpolarization (Lodish et al., 2000; Wright, 2004).

3.1.3. Simulation 3: Influential Factors and Perturbation of Resting V_mem

Simulation 3 explored factors influencing resting V_mem in isolated cells, and also demonstrated the ability for cell V_mem to return to its resting value after a perturbation (Figure 6). As factors, such as membrane permeability to specific ions, ion pump rates, and the influence of environmental ion concentrations, such as high extracellular K⁺ levels, are well known to affect individual cell V_mem (Lodish et al., 2000; Wright, 2004), this simulation (Figure 6) is also a model validation. The simulation was performed on the same cluster used in Simulation 2 (see Figure 5A), with membrane manipulations applied to, and V_mem monitored in, a profile B cell of the cluster. Initial conditions for Simulation 3 were those of the final Simulation 2, with extra/intracellular ion concentrations and V_mem, as listed for the 20 min time point in Table 3.

FIGURE 6

Figure 6. V_mem changes with perturbation to membrane permeability or external ion concentrations. (A) shows membrane permeability to specific ions, which was altered via various forced changes during the simulation. Intervention (a) increased membrane permeability to Na⁺ by a factor of 25 from 1 to 3 s. Intervention (b) increased permeability to K⁺ by a factor of 10 from 13 to 15 s. Intervention (c) increased permeability to Cl⁻ by a factor of 25 from 25 to 27 s, and intervention (d) increased permeability to Ca²⁺ by 50 from 37 to 39 s. In between membrane permeability perturbations, permeability returned to original values. Intervention (e) blocked activity of the Na/K-ATPase pump from 49 to 69 s. Intervention (f) temporarily increased the extracellular concentration of K⁺ by introducing 35 mmol/L KCl at the global boundaries from 79 to 89 s, returning to 5 mmol/L at the boundary after 26 s. (B) shows BESTE-calculated V_mem in comparison to Goldman-derived V_mem for a single cell undergoing the applied interventions, for the situation where instantaneous mixing was assumed in the environmental space (i.e., “cytosol only” simulation). Dashed black line in (B) indicates the value for resting V_mem (−63.7 mV). (C) shows BETSE-calculated V_mem in comparison to Goldman-derived V_mem for a single cell undergoing the same various applied interventions, for the situation where individual extracellular spaces and environmental electrodiffusion were modeled (i.e., “environment modeling” simulation). Dashed black line in (C) indicates the value for resting V_mem (−63.7 mV).

The Goldman equation [equation (22)] was used with membrane permeability and ion concentration values available at each time step to predict V_mem using conventional measures and provide an indicator of expected V_mem.

“Cytosol only” and “environment modeling” are two simulation modes available in BETSE. In “cytosol only” mode, a simple simulation is performed, which assumes instantaneous mixing of fluxes into the environment, where extracellular spaces and free diffusion in the environment are not modeled, and V_mem is calculated assuming the cell is a simple capacitor via the charge inside the cell and the relation $V_{\mathit{mem}}=\frac{1}{\mathit{Cm}}{Q}_{\mathit{cell}}$. The “environment modeling” mode calculates a full extracellular environment using the Maxwell Capacitance Matrix method to solve for voltages, as outlined in the Methods section. The V_mem data presented in Figure 6B is from the “cytosol only” simulation mode, while that from 6C is from the “environment modeling” simulation mode. These two types of simulation modes were compared to illustrate the effect of including extracellular matrix and environmental transport in the bioelectrical model.

Various membrane permeability manipulations, in addition to a block of the Na/K-ATPase pump, and an increase in extracellular K⁺ levels, were simulated (Figure 6). Membrane permeability manipulations effectively simulate the transient opening of an ion channel. The first intervention temporarily increased (from 1 to 3 s) the cell’s membrane permeability to Na⁺ by a factor of 25, leading to a characteristic, pronounced depolarization. The next intervention temporarily increased (from 13 to 15 s) the cell’s membrane permeability to K⁺ by a factor of 10, and generated a characteristic hyperpolarization of V_mem. Subsequent interventions increased the cell’s membrane permeability to Cl⁻ by a factor of 25 from 25 to 27 s, creating an expected depolarization, and increased the membrane permeability to Ca²⁺ by 50 from 37 to 39 s with the expected depolarization effect. The Na/K-ATPase pump activity was blocked from 49 to 69 s, during which time the V_mem for both simulation modes converged at precisely the Goldman V_mem prediction (Figure 6). This result is expected as the Na/K-ATPase pump activity generates an electrogenic current that is not considered in the Goldman analysis (the Goldman equation requires zero net current across the membrane, as discussed previously). Finally, extracellular K⁺ concentration was increased by introducing 35 mmol/L KCl at the global boundaries from 79 to 89 s, which returned to 5 mmol/L at the boundary after the perturbation interval, and resulted in a characteristic, well-known V_mem depolarization (Delamere and Duncan, 1977).

As can be seen in Figure 6, cell V_mem naturally returns to its resting value of −63.7 mV after each perturbation is complete. While overall, V_mem responses for both the “cytosol only” and “environmental modeling” simulation modes captured all main features predicted by the Goldman equation, the “environmental modeling” mode, which includes simulation of extracellular spaces and transport in the environment, showed slower and more complex responses than the “cytosol only” mode.

These results illustrate how physiological circuits implement stability with respect to bioelectric state, as, for example, observed in applications of optogenetics to developing systems (Adams et al., 2013).

3.1.4. Simulation 4: Resting V_mem and Cell Excitability

Simulation set 4 validated the expected function of voltage-gated Na⁺ and K⁺ channels, highlighted the ability for resting V_mem to control cell excitability, and examined the possibility of low voltage-gated K + expression in relation to voltage-gated Na + expression to effect resting V_mem. Simulations were performed on a cluster of 35 cells, which were connected by non-voltage sensitive GJ, and were without TJ. An initialization simulation without voltage-gated channels was run on each cluster to bring cells to equilibrium resting state. Each simulation shown in rows A–C of Figure 7 features cells with different resting V_mem, which is accomplished by altering levels of K⁺ leak channels (altering non-dynamic membrane permeability to K⁺). In simulations A and B of Figure 7, all cells have identical expressions of NaV and KV1.2 channels, with a net maximum membrane permeability of 2667 nm/s and 667 nm/s for NaV and KV1.2 channels, respectively. The resting V_mem of cells in simulation A was −70 mV, while those of simulation B weremuch higher at −18 mV. Simulation C of Figure 7, studied cells with a resting potential of −57 mV and expressions of NaV and KV1.2 channels (homogeneous expressions in the cell population), with a net maximum membrane permeability of 2667 nm/s and 67 nm/s for NaV and KV1.2 channels, respectively. This simulated a deficiency of voltage-gated potassium channels – a phenotype occurring in certain metastatic cancers (Djamgoz, 2014). For each simulation, a forced depolarization is applied to one randomly selected cell of the cluster from a time of 1–200 ms to induce excitable activity.

FIGURE 7

Figure 7. Influence of resting V_mem on cell excitability. Each simulation shown in rows (A–C) features cells with different resting V_mem (as indicated in leftmost panels) which is induced by altering cell levels of K⁺ leak channels (altering non-dynamic membrane permeability to K⁺). In each simulation, all cells have identical NaV and KV1.2 channels, and a forced depolarization is applied to one randomly selected cell of the cluster from a time of 1–200 ms. For cells with the lowest resting potential of −70 mV [row (A)], the forced depolarization leads to the firing of four action-potential-like signals, with excitable activity ceasing with the forced depolarization after 200 ms. For the cluster with the low resting potential of −18 mV [row (B)], the forced depolarization leads to a periodic self-excitation with a period of about 100 ms, which lasts long after the forced depolarization ceases. (C) shows that combined activity of dynamic channels with variable expression levels can generate resting V_mem bi-stability, as for cells with an original resting potential of −57 mV, expression of NaV channels with low expression of KV1.2 transitions the system into a depolarized V_mem of approximately −14 mV after a single forced depolarization.

For cells with the lowest resting potential of −70 mV (Figure 7A), the forced depolarization leads to the firing of four action-potential-like signals, with excitable activity ceasing with the forced depolarization after 200 ms. However, for the cluster with the low resting potential of −18 mV (Figure 7B), the forced depolarization leads to a periodic self-excitation with a period of about 100 ms, which lasts long after the forced depolarization ceases. This demonstrates the well-known expected behavior of cells with resting V_mem higher than the activation threshold of NaV, such as the pacemaker cells of the heart, to enter periodic self-excitations for an indefinite period of time (Roberts and Stirling, 1971; Onganer et al., 2005; Matthews, 2013a). For hyperpolarized cells with resting V_mem of −57 mV with expression of NaV and simulated deficiency of voltage-gated potassium channels (Figure 7B), our simulations indicate that the forced oscillation creates a single action potential, with the resting potential being altered in the long term to a much more depolarized value of −14 mV. These simulations demonstrate both the importance of resting potential in controlling cell excitability, with more depolarized cells showing capacity for self-excitation (Figures 7A,B), and also the capacity for irregular expression of excitable channels to potentially alter the resting V_mem (Figure 7C), which may assist in explaining the sustained depolarization of some cancer cells (Djamgoz, 2014).

3.2. Impacts of Multicellular V_mem Gradients

3.2.1. Simulation 5: Effect of Heterogeneous V_mem on Physiological Properties

Simulation 5 investigated the physiological impacts of a heterogeneous V_mem pattern in a cellular collective. A cluster of 794 cells with a diameter of 375 μm, boundary TJ with a diffusion scaling of 1.0 × 10⁻⁵, and GJ connectivity with a value of ${\mathrm{\beta }}_{\mathit{GJ}}^o$ = 1 × 10⁻⁷was utilized. Initial values for concentrations and voltages in the simulations were those of the final simulation for profile B cells, Table 3. Membrane permeability of cells varied over space in the same pattern and using the same three profiles defined in Figure 5A. The result was a stable pattern of resting V_mem featuring a depolarized spot of cells in the lower left side (V_mem ~ −20 mV) and a hyperpolarized spot of cells in the upper right side of a circular cell cluster (V_mem ~ −60 mV). Each region was surrounded by cells with a mean V_mem of approximately −45 mV (Figure 8A).

FIGURE 8

Figure 8. Differences in V_mem between different cells in a cluster have various biophysical influences on the cluster as a whole. (A) V_mem pattern featuring a depolarized spot of cells in the lower left side and a hyperpolarized spot of cells in the upper right side of a circular cell cluster (A) results in differences in cytosolic Ca²⁺ levels (B), osmotic pressure (C), patterns of voltage-sensitive gap junction conductivity [(D), where black is open and white is closed], and small interior environmental voltage gradients (E). Patterns of V_mem also influence the cytosolic and extracellular concentration of a negatively charged anionic compound (F), induce strong nano-scale electric fields between gap junctions (G), and generate a characteristic long-range pattern of total ionic current density and macroscopic electric field (H).

The presence of regional V_mem differences was found to have various influences on the cluster as a whole. Heterogeneous V_mem induced differences in cytosolic Ca²⁺ levels (Figure 8B) in a manner inversely proportional to cell V_mem, with the most hyperpolarized cells having cytosolic Ca²⁺ of over 150 nmol/L while the most depolarized contained <60 nmol/L. By contrast, a hypothetical negatively charged anionic signaling molecule develops a cytosolic concentration profile in direct correspondence to V_mem values (inverse to that of cationic Ca²⁺), but due to the presence of TJ, which enable the formation of extracellular voltages due to charge internalization (Figure 8F), the anionic substance concentrates in extracellular spaces around hyperpolarized cells (Figure 8F).

Heterogeneous V_mem was also seen to produce significant osmotic pressure differences between cells of different resting potential, with more hyperpolarized V_mem leading to lower osmotic pressure than more depolarized V_mem (Figure 8C). This is consistent with expectations, as in simulation 5 more hyperpolarized cells have higher K⁺ leak channels, which means more K⁺ is moving out of the cell and into extracellular spaces with expected water movement from the cytosol to the extracellular space to compensate for movement of salt (i.e., lowering of osmotic pressure). By contrast, depolarized cells of the simulation have higher levels of Na⁺ leak permeability, which means more Na⁺ is moving from the extracellular space to the cytosol with expected water movement from the extracellular space to the cell to compensate (increase of intracellular osmotic pressure). Depending on the mechanical properties of cells, these osmotic pressures may translate into cell volume changes and hydrostatic pressures and pressure gradients (body forces).

Voltage-sensitive Gap junctions connecting cells responded to voltage gradients created by V_mem, closing to minimum conductivity values and isolating the two regions of differential V_mem from the remainder of the cluster (Figure 8D). The V_mem pattern in this example generated electric fields of up to ~6.5 × 10⁵ V/m acting (over short spatial distances of 26 nm) between interfacing cell membranes of GJ networked cells (Figure 8G).

Heterogeneous V_mem in a GJ networked multicellular cluster also was found to induce a long-range pattern of total ionic current density up to 60 μA/cm² in magnitude (Figure 8H).

We conclude that stable patterns of resting V_mem have numerous, significant impacts on the cluster as a whole, altering concentration profiles of key signaling moieties, inducing physiological pressures and forces, and establishing long-range patterns of ion transport and macroscopic electric field.

3.3. Factors Influencing Resting V_mem in Networked Cells

3.3.1. Simulation 6: Gap Junction Connectivity

To understand group dynamics and the dynamics of bioelectric states in an electrically coupled tissue, this simulation explored the influence of GJ connectivity on cell resting V_mem for a small group of cells (encircled in Figure 9) with 15× increased Na⁺ membrane permeability (simulating an increased expression of an open Na⁺ ion channel). The multicellular cluster contained 794 cells, and had a diameter of 375 μm. Cells were connected by GJ at interfacing membranes. The simulation began with values for intra/extracellular concentrations and V_mem obtained after a 20 minute initialization simulation, which were similar to those listed in Table 3’s 20 min time point for profile B cells. All cells began with identical membrane permeability profiles with values of profile B cells as listed in Figure 5A.

FIGURE 9

Figure 9. Gap junction connectivity has a strong affect on the ability for cells with altered properties to manifest different bioelectrical states in the collective. A small patch of cells developing 15× increased Na⁺ permeability transitions those cells into a new resting V_mem for cells with low GJ connectivity (A) but not for cells with high GJ connectivity (B). (C) shows the time dependence of Na membrane permeability for the affected patch of cells, while (D) shows the time course of V_mem for a cell in the affected patch of cells for high and low GJ connectivity clusters.

Cells in a first simulation (low GJ connectivity) had an intercellular (GJ mediated) free-diffusion scaling of ${\mathrm{\beta }}_{\mathit{GJ}}^o$ = 1.0 × 10⁻⁷. For a cell with 1.0 × 10⁵ GJ in total, this corresponds to an individual GJ conductivity of approximately 68 pS.

Cells in a second simulation (low GJ connectivity) had an intercellular (GJ mediated) diffusion scaling constant of ${\mathrm{\beta }}_{\mathit{GJ}}^o$ = 1.0 × 10⁻⁶. For a cell with 1.0 × 10⁵ GJ in total, this corresponds to an individual GJ conductivity of approximately 6.8 pS.

For both high and low GJ connectivity simulations, at t = 1.0 s of the simulation Na⁺ membrane permeability of a small patch of cells (circled in Figures 9A,B) was increased by 15× and remained increased for the duration of the simulation (Figure 9C). This simulates the increased expression, or opening, of a Na⁺ ion channel in this small patch of cells, but not in the remaining cells of the cluster.

Effects on V_mem vary significantly between cells with high or low GJ connectivity (Figure 9). For a cluster with low GJ connectivity, the 15× increase in Na⁺ permeability leads to approximately 40 mV depolarization of V_mem, which remains stable as a new resting V_mem state divergent from that of surrounding cells (Figure 9D). However, the cluster with high GJ connectivity shows only 10 mV depolarization in V_mem with the 15× increase in Na⁺ permeability (Figure 9D).

We conclude that the characteristics of GJ connectivity in a cluster have a significant influence in specifying resting V_mem for cells with heterogeneous ion channel characteristics, and may, therefore, be expected to play important roles in morphogenesis and the development of cancer, which both require the development of differential V_mem states from a homogeneous collective.

3.3.2. Simulation 7: Tight Junction Connectivity

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⁻⁵ 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).

FIGURE 10

Figure 10. A trans-boundary voltage difference (trans-epithelial potential; TEP) appears across the outer boundary of cell clusters when ion transport between extracellular spaces at the cluster boundary is limited by simulated tight junctions (B,C). This trans-boundary voltage does not occur when ion transport in extracellular spaces is similar to free-diffusion values for ions (A). With TJ present, wounding leads to a characteristic bioelectrical cluster state (C) where ion current flows out of the wounded area, returning to the cluster at both sides of the wound, and establishing a long-range (100 s of micrometers) current flow and macroscopic electric field in the cluster (C).

FIGURE 11

Figure 11. The TEP appearing across the outer boundary of cell clusters, when ion transport between extracellular spaces at the cluster boundary is limited by simulated tight junctions, was indicated to form because of charge and voltage build up in the extracellular environment (A). The TEP strongly polarizes individual cells in the outer layer of the cluster (B), and results in an electric field at the cluster boundary. The strong voltage gradients of the TEP were predicted to create spatial patterns of GJ voltage gating, whereby cells in the outer layer and cluster interior form two networked groups, but cells between the outer layer and cluster interior are separated by closed and minimally conductive GJ [(C), black represents fully open GJ, white represents fully closed GJ]. Patterns of global ionic current density (not shown) and electroosmotic fluid flow (E) are strongest in the cluster across the TEP. The strong electric field and electroosmotic flows at the cluster boundary are predicated to result in lateral self-electrophoresis/electroosmosis of ion pumps and channels in the membranes (D), which further alters current, flow, TEP, and V_mem patterns of the cluster.

The mere presence of TJ, even in the absence of inhomogeneous membrane channel distribution, current, or electroosmotic flow, was seen to alter cell V_mem 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 V_mem 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.

3.3.3. Simulation 8: Spontaneous V_mem Patterning

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 V_mem 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 V_mem sensitivity of ion channels and GJ, the non-linear relationship between channel/GJ activity and V_mem, 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 V_mem gradient in a cluster of voltage-sensitive GJ networked cells with TJ (Figure 12). In its initial state, the cluster had a small V_mem 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 V_mem gradient (Figure 12).

FIGURE 12

Figure 12. Example feedback mechanism generating V_mem patterning in a cluster of voltage-sensitive, GJ networked cells with TJ, which occurs due to the presence of a charged channel gating ligand. In its initial state, a cell cluster has a small V_mem asymmetry due to increase in K⁺ leak channels for a small cluster of cells in the upper right side, leading to a minor resting V_mem difference of approximately 10 mV [Initial State (A)]. This initial state also corresponds to small related asymmetries in environmental voltage [Initial State (B)], anionic ligand concentration [Initial State (C)], and no changes to GJ open/closed state [Initial State (D), where blue = open, white = closed]. However, the presence of an electrodiffusing cationic gating ligand capable of opening a sodium channel based on its extracellular concentration (analogous to acetylcholine), in combination with TJ and GJ activity, generates a strong axial V_mem polarity on account of the positive feedback mechanism outlined in (E).

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 V_mem 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 V_mem becoming more depolarized. As voltage gradients develop between cells, GJ close, further reinforcing the development of two regions of highly distinguished V_mem (Figure 12B).

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

4. Discussion

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 V_mem on the basis of simulated electrodiffusing “reporter dye” intra- and extracellular concentrations, showed remarkable correspondence to BETSE direct-calculated V_mem. Likewise, comparison between BETSE direct-calculated and Goldman-derived V_mem 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/cm² were within the range of typically observed experimental currents 1–500 μA/cm², 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 V_mem 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 V_mem 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 V_mem 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 V_mem 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 V_mem from a hyperpolarized (−57 mV) to depolarized resting V_mem state (−14 mV) after a single transient depolarization event, which is consistent with the abnormal resting V_mem observed in some cancers (−20 to −5 mV) (Binggeli and Weinstein, 1986; Djamgoz, 2014). A similar form of resting V_mem 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 V_mem can exert. One-way V_mem can exert an influence on downstream patterning mechanisms is by altering the spatial distribution of important chemical signaling molecules (such as Ca²⁺, serotonin, or glutamate) in relation to cell V_mem, 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 V_mem 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 V_mem differences can passively induce differences in cytosolic Ca²⁺ levels (Figure 8B) in a manner inversely related to cell V_mem. As Ca²⁺ is an important secondary messenger, and intracellular Ca²⁺ levels are involved in calcium-induced-calcium-release (CICR) and ion channel gating (e.g., Ca²⁺ gated K⁺ channels), these differences in the cytosolic Ca²⁺ concentration with V_mem may produce significant downstream effects on cell state (including enzyme function, gene expression, and apoptosis) and further evolution of bioelectric pattern (via effects on Ca²⁺ gated ion channels). Similarly, effects on the concentration profile of a negatively charged signaling molecule were in direct correspondence to V_mem values (inverse to that of cationic Ca²⁺), which further illustrates how cell state can be influenced in divergent ways as by simply exhibiting different values of resting V_mem signaling molecules with different charge have opposite changes to their concentration profiles. Further V_mem 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 V_mem states (McCaig et al., 2005; Tseng and Levin, 2013; Levin, 2014), clear comprehension of the factors involved in specifying cell resting V_mem 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 V_mem 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 V_mem 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 V_mem 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 V_mem 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 V_mem states in the collective (Figure 9). Since research indicates that resting V_mem 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 V_mem 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 V_mem pattern.

The involvement of GJ- and TJ-related factors in the specification of multicellular V_mem 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 V_mem and asymmetry in developing clusters (Figures 9 and 10). In cancer, which is characterized by cells with depolarized V_mem, loss of GJ communication would allow cells with different membrane channel states to express dramatically altered V_mem 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 V_mem, which in turn can affect voltage-sensitive channels and electrical synapses. Furthermore, the presence of ionic messengers, including Ca²⁺ or anionic serotonin or glutamate, are charged molecules subject to movement in electrochemical gradients, but can also alter V_mem 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 V_mem 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 V_mem 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).

Author Contributions

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.

Conflict of Interest Statement

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.

Acknowledgments

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.

Funding

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).

Supplementary Material

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

References

Adamatzky, A., Holley, J., Dittrich, P., Gorecki, J., De Lacy Costello, B., Zauner, K. P., et al. (2012). On architectures of circuits implemented in simulated belousov-zhabotinsky droplets. BioSystems 109, 72–77. doi: 10.1016/j.biosystems.2011.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Adamatzky, A., and Jones, J. (2011). On electrical correlates of physarum polycephalum spatial activity: can we see physarum machine in the dark? Bioph. Rev. Lett. 6, 29–57. doi:10.1142/S1793048011001257

CrossRef Full Text | Google Scholar

Adams, D., Tseng, A., and Levin, M. (2013). Light-activation of the archaerhodopsin H(+)-pump reverses age-dependent loss of vertebrate regeneration: sparking system-level controls in vivo. Biol. Open 2, 306–313. doi:10.1242/bio.20133665

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, D. S., Masi, A., and Levin, M. (2007). H+ pump-dependent changes in membrane voltage are an early mechanism necessary and sufficient to induce Xenopus tail regeneration. Development 134, 1323–1335. doi:10.1242/dev.02812

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, D. S., Uzel, S. G., Akagi, J., Wlodkowic, D., Andreeva, V., Yelick, P. C., et al. (2016). Bioelectric signalling via potassium channels: a mechanism for craniofacial dysmorphogenesis in kcnj2-associated Andersen-Tawil syndrome. J. Physiol. 594, 3245–3270. doi:10.1113/JP271930

PubMed Abstract | CrossRef Full Text | Google Scholar

Altizer, A., Moriarty, L., Bell, S., Schreiner, C., Scott, W., and Borgens, R. (2001). Endogenous electric current is associated with normal development of the vertebrate limb. Dev. Dyn. 221, 391–401. doi:10.1002/dvdy.1158

PubMed Abstract | CrossRef Full Text | Google Scholar

Andreev, V. (2013). Cytoplasmic electric fields and electroosmosis: possible solution for the paradoxes of the intracellular transport of biomolecules. PLoS ONE 8:e61884. doi:10.1371/journal.pone.0061884

PubMed Abstract | CrossRef Full Text | Google Scholar

Arcangeli, A., Crociani, O., Lastraioli, E., Masi, A., Pillozzi, S., and Becchetti, A. (2009). Targeting ion channels in cancer: a novel frontier in antineoplastic therapy. Curr. Med. Chem. 16, 66–93. doi:10.2174/092986709787002835

PubMed Abstract | CrossRef Full Text | Google Scholar

Barel, O., Shalev, S. A., Ofir, R., Cohen, A., Zlotogora, J., Shorer, Z., et al. (2008). Maternally inherited birk barel mental retardation dysmorphism syndrome caused by a mutation in the genomically imprinted potassium channel KCNK9. Am. J. Hum. Genet. 83, 193–199. doi:10.1016/j.ajhg.2008.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazant, M., and Squires, T. (2004). Induced-charge electrokinetic phenomena: theory and microfluidic applications. Phys. Rev. Lett. 92, 066101. doi:10.1103/PhysRevLett.92.066101

PubMed Abstract | CrossRef Full Text | Google Scholar

Beane, W., Morokuma, J., Lemire, J., and Levin, M. (2013). Bioelectric signaling regulates head and organ size during planarian regeneration. Development 140, 313–322. doi:10.1242/dev.086900

PubMed Abstract | CrossRef Full Text | Google Scholar

Beane, W. S., Morokuma, J., Adams, D. S., and Levin, M. (2011). A chemical genetics approach reveals H,K-ATPase-mediated membrane voltage is required for planarian head regeneration. Chem. Biol. 18, 77–89. doi:10.1016/j.chembiol.2010.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Beard, D., and Qian, H. (2007). Relarelation between thermodynamic driving force and one-way fluxes in reversible processes. PLoS ONE 2:e144. doi:10.1371/journal.pone.0000144

CrossRef Full Text | Google Scholar

Bendahhou, S., Donaldson, M. R., Plaster, N. M., Tristani-Firouzi, M., Fu, Y. H., and Ptacek, L. J. (2003). Defective potassium channel Kir2.1 trafficking underlies Andersen-Tawil syndrome. J. Biol. Chem. 278, 51779–51785. doi:10.1074/jbc.M310278200

PubMed Abstract | CrossRef Full Text | Google Scholar

Berridge, M. (2014). Ion channels. Cell Signal. Biol. 3, 1–72.

Google Scholar

Binggeli, R., and Weinstein, R. (1986). Membrane potentials and sodium channels: hypotheses for growth regulation and cancer formation based on changes in sodium channels and gap junctions. J. Theor. Biol. 123, 377–401. doi:10.1016/S0022-5193(86)80209-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Blackiston, D., Adams, D. S., Lemire, J. M., Lobikin, M., and Levin, M. (2011). Transmembrane potential of glycl-expressing instructor cells induces a neoplastic-like conversion of melanocytes via a serotonergic pathway. Dis. Model Mech. 4, 67–85. doi:10.1242/dmm.005561

PubMed Abstract | CrossRef Full Text | Google Scholar

Blackiston, D., McLaughlin, K., and Levin, M. (2009). Bioelectric controls of cell proliferation: ion channels, membrane voltage and the cell cycle. Cell Cycle 8, 3527–3536. doi:10.4161/cc.8.21.9888

PubMed Abstract | CrossRef Full Text | Google Scholar

Borgens, R., Callahan, L., and Rouleau, M. (1987). Anatomy of axolotl flank integument during limb bud development with special reference to a transcutaneous current predicting limb formation. J. Exp. Zool. 244, 203–214. doi:10.1002/jez.1402440204

PubMed Abstract | CrossRef Full Text | Google Scholar

Borgens, R. B. (1984). Are limb development and limb regeneration both initiated by an integumentary wounding? A hypothesis. Differentiation 28, 87–93. doi:10.1111/j.1432-0436.1984.tb00270.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Boutillier, A. L., Kienlen-Campard, P., and Loeffler, J. P. (1999). Depolarization regulates cyclin d1 degradation and neuronal apoptosis: a hypothesis about the role of the ubiquitin/proteasome signalling pathway. Eur. J. Neurosci. 11, 441–448. doi:10.1046/j.1460-9568.1999.00451.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bower, J. M., and Beeman, D. (2007). Constructing realistic neural simulations with genesis. Methods Mol. Biol. 401, 103–125. doi:10.1007/978-1-59745-520-6_7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowman, C., and Baglioni, A. (1984). Application of the Goldman-Hodgkin-Katz current equation to membrane current-voltage data. J. Theor. Biol. 108, 1–29. doi:10.1016/S0022-5193(84)80165-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bronstein-Sitton, N. (2004). Regulating the immune response: the unexpected role of ion channels. Modulator 18, 1–5.

Google Scholar

Cervera, J., Alcaraz, A., and Mafe, S. (2014). Membrane potential bistability in nonexcitable cells as described by inward and outward voltage-gated ion channels. J. Phys. Chem. B 118, 12444–12450. doi:10.1021/jp508304h

PubMed Abstract | CrossRef Full Text | Google Scholar

Cervera, J., Alcaraz, A., and Mafe, S. (2016). Bioelectrical signals and ion channels in the modeling of multicellular patterns and cancer biophysics. Sci. Rep. 6, 20403. doi:10.1038/srep20403

PubMed Abstract | CrossRef Full Text | Google Scholar

Cervera, J., Manzanares, J. A., and Mafe, S. (2015). Electrical coupling in ensembles of nonexcitable cells: modeling the spatial map of single cell potentials. J. Phys. Chem. B 119, 2968–2978. doi:10.1021/jp512900x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chernet, B., and Levin, M. (2013). Transmembrane voltage potential is an essential cellular parameter for the detection and control of tumor development in a Xenopus model. Dis. Model. Mech. 6, 595–607. doi:10.1242/dmm.010835

PubMed Abstract | CrossRef Full Text | Google Scholar

Chernet, B., and Levin, M. (2014). Transmembrane voltage potential of somatic cells controls oncogene-mediated tumorigenesis at long-range. Oncotarget 5, 3287–3306. doi:10.18632/oncotarget.1935

PubMed Abstract | CrossRef Full Text | Google Scholar

Chernet, B. T., Fields, C., and Levin, M. (2015). Long-range gap junctional signaling controls oncogene-mediated tumorigenesis in Xenopus laevis embryos. Front. Physiol. 5:519. doi:10.3389/fphys.2014.00519

PubMed Abstract | CrossRef Full Text | Google Scholar

Clements, J., Clayton, R., and Arlon, T. (1975). Computation of the capacitance matrix for systems of dielectric-coated cylindrical conductors. IEEE Trans. Electromagn. Compat. ECM-17, 238–249. doi:10.1109/TEMC.1975.303430

CrossRef Full Text | Google Scholar

Costa, F., Emilio, M., Fernandes, P., Gil Ferreira, H., and Ferreira, K. G. (1989). Determination of ionic permeability coefficients of the plasma membrane of Xenopus laevis oocytes under voltage clamp. J. Physiol. 413, 199–211. doi:10.1113/jphysiol.1989.sp017649

PubMed Abstract | CrossRef Full Text | Google Scholar

Cross, M., and Hohenberg, P. (1993). Pattern formation outside of equilibrium. Rev. Mod. Phys. 65, 851. doi:10.1103/RevModPhys.65.851

CrossRef Full Text | Google Scholar

De Loof, A. (1985). The cell as a miniature electrophoresis chamber. Comp. Biochem. Physiol. 80A, 453–459. doi:10.1016/0300-9629(85)90397-4

CrossRef Full Text | Google Scholar

Delamere, N., and Duncan, G. (1977). A comparison of ion concentrations, potentials and conductances of amphibian, bovine and cephalopod lenses. J. Physiol. 272, 167–186. doi:10.1113/jphysiol.1977.sp012039

PubMed Abstract | CrossRef Full Text | Google Scholar

Djamgoz, M. (2014). Biophysics of cancer: cellular excitability (“celex”) hypothesis of metastasis. J. Clin. Exp. Oncol. S1, 005. doi:10.4172/2324-9110.S1-005

CrossRef Full Text | Google Scholar

Doursat, R., and Sanchez, C. (2014). Growing fine-grained multicellular robots. Soft Rob. 1, 110–121. doi:10.1089/soro.2014.0014

CrossRef Full Text | Google Scholar

Emmons-Bell, M., Durant, F., Hammelman, J., Bessonov, N., Volpert, V., Morokuma, J., et al. (2015). Gap junctional blockade stochastically induces different species-specific head anatomies in genetically wild-type girardia dorotocephala flatworms. Int. J. Mol. Sci. 16, 27865–27896. doi:10.3390/ijms161126065

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, J., Sun, X., Moore, L., White, T., Brink, P., and Mathia, R. (2011). Lens intracellular hydrostatic pressure is generated by the circulation of sodium and modulated by gap junction coupling. J. Gen. Physiol. 137, 507–520. doi:10.1085/jgp.201010538

PubMed Abstract | CrossRef Full Text | Google Scholar

Geard, N., and Willadsen, K. (2009). Dynamical approaches to modeling developmental gene regulatory networks. Birth Defects Res. C Embryo Today 87, 131–142. doi:10.1002/bdrc.20150

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodenough, D., and Paul, D. (2009). Gap junctions. Cold Spring Harb. Perspect. Biol. 1, a002576. doi:10.1101/cshperspect.a002576

PubMed Abstract | CrossRef Full Text | Google Scholar

Gow, N. A., and Morris, B. M. (1995). The electric fungus. Bot. J. Scotl. 47, 263–277. doi:10.1080/03746609508684833

CrossRef Full Text | Google Scholar

Hamill, O., Huguenard, J., and Prince, D. (1991). Patch-clamp studies of voltage-gated currents in identified neurons of the rat cerebral cortex. Cereb. Cortex 1, 48–61. doi:10.1093/cercor/1.1.48

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, A., Spray, D., and Bennett, M. (1983). Control of intercellular communication by voltage dependence of gap junctional conductance. J. Neurosci. 3, 79–100.

PubMed Abstract | Google Scholar

Hay, J., and Geddes, D. (1985). Transepithelial potential difference in cystic fibrosis. Thorax 40, 493–496. doi:10.1136/thx.40.7.493

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinzel, T. (2008). “Chapter C: capacitance matrix and electrostatic energy,” in Mesoscopic Electronics in Solid State Nanostructures, Second Edition (Weinheim: Wiley), 349–352.

Google Scholar

Hotary, K., and Robertson, F. (1994). Endogenous electric currents and voltage gradients in Xenopus embryos and the consequences of their disruption. Dev. Biol. 166, 789–800. doi:10.1006/dbio.1994.1357

CrossRef Full Text | Google Scholar

Hotary, K., and Robinson, K. (1990). Endogenous electrical currents and the resultant voltage gradients in the chick embryo. Dev. Biol. 140, 149–160. doi:10.1016/0012-1606(90)90062-N

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaffe, L. (1981). The role of ionic currents in establishing developmental pattern. Philos. Trans. R. Soc. Lond. B Biol. Sci. 295, 553–566. doi:10.1098/rstb.1981.0160

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamm, R. D., and Bashir, R. (2014). Creating living cellular machines. Ann. Biomed. Eng. 42, 445–459. doi:10.1007/s10439-013-0902-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Law, R., and Levin, M. (2015). Bioelectric memory: modeling resting potential bistability in amphibian embryos and mammalian cells. Theor. Biol. Med. Model. 12, 22. doi:10.1186/s12976-015-0019-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Leithe, E., Sirnes, S., Omori, Y., and Rivedal, E. (2006). Downregulation of gap junctions in cancer cells. Crit. Rev. Oncog. 12, 225–256. doi:10.1615/CritRevOncog.v12.i3-4.30

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, M. (2012). Morphogenetic fields in embryogenesis, regeneration, and cancer: non-local control of complex patterning. BioSystems 109, 243–261. doi:10.1016/j.biosystems.2012.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, M. (2014). Molecular bioelectricity: how endogenous voltage potentials control cell behavior and instruct pattern regulation in vivo. Mol. Biol. Cell 25, 3835–3850. doi:10.1091/mbc.E13-12-0708

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, M., Buznikov, G., and Lauder, J. (2006). Of minds and embryos: left-right asymmetry and the serotonergic controls of pre-neural morphogenesis. Dev. Neurosci. 28, 171–185. doi:10.1159/000091915

PubMed Abstract | CrossRef Full Text | Google Scholar

Levin, M., and Stephenson, C. (2012). Regulation of cell behavior and tissue patterning by bioelectrical signals: challenges and opportunities for biomedical engineering. Annu. Rev. Biomed. Eng. 14, 295–323. doi:10.1146/annurev-bioeng-071811-150114

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobikin, M., Chernet, B., Lobo, D., and Levin, M. (2012). Resting potential, oncogene-induced tumorigenesis, and metastasis: the bioelectric basis of cancer in vivo. Phys. Biol. 9, 065002. doi:10.1088/1478-3975/9/6/065002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobo, D., and Levin, M. (2015). Inferring regulatory networks from experimental morphological phenotypes: a computational method reverse-engineers planarian regeneration. PLoS Comput. Biol. 4:e1004295. doi:10.1371/journal.pcbi.1004295

PubMed Abstract | CrossRef Full Text | Google Scholar

Lodish, H., Berk, A., Zipursky, S., Matsudaira, P., Baltimore, D., and Darnell, J. (2000). “Section 15.4: intracellular ion environment and membrane electric potential,” in Molecular Cell Biology 4th Edition (New York: W.H. Freeman and Co).

Google Scholar

Marsh, G., and Beams, H. W. (1952). Electrical control of morphogenesis in regenerating Dugesia tigrina. I. Relation of axial polarity to field strength. J. Cell. Comp. Physiol. 39, 191–213. doi:10.1002/jcp.1030390203

CrossRef Full Text | Google Scholar

Matthews, G. (2013a). “Appendix A: derivation of the Nernst equation,” in Cellular Physiology of Nerve and Muscle, Fourth Edition (Berlin: Wiley), 208–211.

Google Scholar

Matthews, G. (2013b). “Appendix B: derivation of the Goldman equation,” in Cellular Physiology of Nerve and Muscle, Fourth Edition (Berlin: Wiley), 212–215.

Google Scholar

Mazzanti, M., Bustamante, J., and Oberleithner, H. (2001). Electrical dimension of the nuclear envelope. Physiol. Rev. 81, 1–19.

PubMed Abstract | Google Scholar

McCaig, C., Rajnicek, A., Song, B., and Zhao, M. (2005). Controlling cell behavior electrically: current views and future potential. Physiol. Rev. 85, 943–978. doi:10.1152/physrev.00020.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

McCaig, C. D. (1990). Nerve branching is induced and oriented by a small applied electric field. J. Cell. Sci. 95(Pt 4), 605–615.

PubMed Abstract | Google Scholar

McLaughlin, S., and Poo, M. (1981). The role of electroosmosis in the electric-field induced movement of charged macromolecules on the surfaces of cells. Biophys. J. 34, 85–93. doi:10.1016/S0006-3495(81)84838-2

CrossRef Full Text | Google Scholar

Millman, J. K., and Aivazis, M. (2011). Python for scientists and engineers. Comput. Sci. Eng. 13, 9–12. doi:10.1109/MCSE.2011.36

CrossRef Full Text | Google Scholar

Milnor, J. (1985). On the concept of attractor. Commun. Math. Phys. 99, 177–195. doi:10.1007/BF01212280

CrossRef Full Text | Google Scholar

Morokuma, J., Blackiston, D., Adams, D. S., Seebohm, G., Trimmer, B., and Levin, M. (2008). Modulation of potassium channel function confers a hyperproliferative invasive phenotype on embryonic stem cells. Proc. Natl. Acad. Sci. U.S.A. 105, 16608–16613. doi:10.1073/pnas.0808328105

PubMed Abstract | CrossRef Full Text | Google Scholar

Munzer, J., Daly, S., Jewell-Motz, E., Lingrel, J., and Blostein, R. (1994). Tissue- and isoform-specific kinetic behavior of the Na,K-ATPase. J. Biol. Chem. 269, 16668–16676.

PubMed Abstract | Google Scholar

Mustard, J., and Levin, M. (2014). Bioelectrical mechanisms for programming growth and form: taming physiological networks for soft body robotics. Soft Rob. 1, 169–191. doi:10.1089/soro.2014.0011

CrossRef Full Text | Google Scholar

Nau, C. (2008). Voltage-gated ion channels. Handb. Exp. Pharmacol. 182, 85–92. doi:10.1007/978-3-540-74806-9_4

CrossRef Full Text | Google Scholar

Nelson, M. (2004). “Electrophysiological models,” in Databasing the Brain: From Data to Knowledge, ed. S. H. Koslow (New York: Wiley), 285–301.

Google Scholar

Ng, S., Chin, C., Lau, Y., Luo, J., Wong, C., Bian, Z., et al. (2010). Role of voltage-gated potassium channels in the fate determination of embryonic stem cells. J. Cell. Physiol. 224, 165–177. doi:10.1002/jcp.22113

PubMed Abstract | CrossRef Full Text | Google Scholar

Nuccitelli, R. (1992). Endogenous ionic currents and DC electric fields in multicellular animal tissues. Bioelectromagnetics S1, 147–157. doi:10.1002/bem.2250130714

PubMed Abstract | CrossRef Full Text | Google Scholar

Nuccitelli, R. (2003a). Endogenous electric fields in embryos during development, regeneration and wound healing. Radiat. Prot. Dosimetry 106, 375–383. doi:10.1093/oxfordjournals.rpd.a006375

CrossRef Full Text | Google Scholar

Nuccitelli, R. (2003b). A role for endogenous electric fields in wound healing. Curr. Top. Dev. Biol. 58, 1–26. doi:10.1016/S0070-2153(03)58001-2

CrossRef Full Text | Google Scholar

Onganer, P., Seckl, M., and Djamgoz, M. (2005). Neuronal characteristics of small-cell lung cancer. Br. J. Cancer 93, 1197–1201. doi:10.1038/sj.bjc.6602857

PubMed Abstract | CrossRef Full Text | Google Scholar

Oviedo, N., Morokuma, J., Walentek, P., Kema, I., Gu, M., Ahn, J., et al. (2010). Long-range neural and gap junction protein-mediated cues control polarity during planarian regeneration. Dev. Biol. 339, 188–199. doi:10.1016/j.ydbio.2009.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Pai, V., Aw, S., Shomrat, T., Lemire, J., and Levin, M. (2012). Transmembrane voltage potential controls embryonic eye patterning in Xenopus laevis. Development 139, 313–323. doi:10.1242/dev.073759

PubMed Abstract | CrossRef Full Text | Google Scholar

Pai, V., Lemire, J., Paré, J., Lin, G., Chen, Y., and Levin, M. (2015). Endogenous gradients of resting potential instructively pattern embryonic neural tissue via notch signaling and regulation of proliferation. J. Neurosci. 35, 4366–4385. doi:10.1523/JNEUROSCI.1877-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Palacios-Prado, N., and Bukauskas, F. F. (2009). Heterotypic gap junction channels as voltage-sensitive valves for intercellular signaling. Proc. Natl. Acad. Sci. U.S.A. 106, 14855–14860. doi:10.1073/pnas.0901923106

PubMed Abstract | CrossRef Full Text | Google Scholar

Pekar, M. (2015). The thermodynamic driving force for kinetics in general and enzyme kinetics in particular. Chemphyschem 16, 884–885. doi:10.1002/cphc.201402778

PubMed Abstract | CrossRef Full Text | Google Scholar

Perathoner, S., Daane, J. M., Henrion, U., Seebohm, G., Higdon, C. W., Johnson, S. L., et al. (2014). Bioelectric signaling regulates size in zebrafish fins. PLoS Genet. 10:e1004080. doi:10.1371/journal.pgen.1004080

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranjan, R., Khazen, G., Gambazzi, L., Ramaswamy, S., Hill, S., Schürmann, F., et al. (2011). Channelpedia: an integrative and interactive database for ion channels. Front. Neuroinform. 5:1–8. doi:10.3389/fninf.2011.00036

PubMed Abstract | CrossRef Full Text | Google Scholar

Raspopovic, J., Marcon, L., Russo, L., and Sharpe, J. (2014). Modeling digits. Digit patterning is controlled by a bmp-sox9-Wnt turing network modulated by morphogen gradients. Science 345, 566–570. doi:10.1126/science.1252960

PubMed Abstract | CrossRef Full Text | Google Scholar

Reingruber, J., and Holcman, D. (2014). Computational and mathematical methods for morphogenetic gradient analysis, boundary formation and axonal targeting. Semin. Cell Dev. Biol. 35, 189–202. doi:10.1016/j.semcdb.2014.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Roberts, A., and Stirling, C. (1971). The properties and propagation of a cardiac-like impulse in the skin of young tadpoles. Z. Vgl. Physiol. 71, 295–310. doi:10.1007/BF00298141

CrossRef Full Text | Google Scholar

Schafer, M. (2006). “Chapter 4: finite volume methods,” in Computational Engineering – Introduction to Numerical Methods (Berlin; Heidelberg: Springer), 77–105.

Google Scholar

Shi, R., and Borgens, R. (1995). Three-dimensional gradients of voltage during development of the nervous system as invisible coordinates for the establishment of embryonic pattern. Dev. Dyn. 202, 101–114. doi:10.1002/aja.1002020202

PubMed Abstract | CrossRef Full Text | Google Scholar

Slack, J. (2014). Establishment of spatial pattern. Wiley Interdiscip. Rev. Dev. Biol. 3, 379–388. doi:10.1002/wdev.144

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, N., and Crampin, E. (2004). Development of models of active ion transport for whole-cell modelling: cardiac sodium-potassium pump as a case study. Prog. Biophys. Mol. Biol. 85, 387–405. doi:10.1016/j.pbiomolbio.2004.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Soler, A., Miller, R., Laughlin, K., NZ, C., Klurfeld, D., and Mullin, J. (1999). Increased tight junctional permeability is associated with the development of colon cancer. Carcinogenesis 20, 1425–1431. doi:10.1093/carcin/20.8.1425

PubMed Abstract | CrossRef Full Text | Google Scholar

Sprunger, L., Stewig, N., and O’Grady, S. (1996). Effects of charybdotoxin on K+ channel (KV1.2) deactivation and inactivation kinetics. Eur. J. Pharmacol. 314, 357–364. doi:10.1016/S0014-2999(96)00556-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundelacruz, S., Levin, M., and Kaplan, D. (2008). Membrane potential controls adipogenic and osteogenic differentiation of mesenchymal stem cells. PLoS ONE 3:e3737. doi:10.1371/journal.pone.0003737

PubMed Abstract | CrossRef Full Text | Google Scholar

Trosko, J. E. (2007). Gap junctional intercellular communication as a biological “rosetta stone” in understanding, in a systems biological manner, stem cell behavior, mechanisms of epigenetic toxicology, chemoprevention and chemotherapy. J. Membr. Biol. 218, 93–100. doi:10.1007/s00232-007-9072-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Tseng, A., and Levin, M. (2013). Cracking the bioelectric code: probing endogenous ionic controls of pattern formation. Commun. Integr. Biol. 6, 1–8. doi:10.4161/cib.22595

PubMed Abstract | CrossRef Full Text | Google Scholar

Tseng, A. S., Beane, W. S., Lemire, J. M., Masi, A., and Levin, M. (2010). Induction of vertebrate regeneration by a transient sodium current. J. Neurosci. 30, 13192–13200. doi:10.1523/JNEUROSCI.3315-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Tseng, A. S., Carneiro, K., Lemire, J. M., and Levin, M. (2011). HDAC activity is required during Xenopus tail regeneration. PLoS ONE 6:e26382. doi:10.1371/journal.pone.0026382

PubMed Abstract | CrossRef Full Text | Google Scholar

Turing, A. (1952). The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B Biol. Sci. 237, 37–72. doi:10.1098/rstb.1952.0012

CrossRef Full Text | Google Scholar

Vandenberg, L., Morrie, R., and Adams, D. (2011). V-ATPase-dependent ectodermal voltage and ph regionalization are required for craniofacial morphogenesis. Dev. Dyn. 240, 1889–1904. doi:10.1002/dvdy.22685

PubMed Abstract | CrossRef Full Text | Google Scholar

Veech, R., Kashiwaya, Y., and King, T. (1995). The resting membrane potential of cells are measures of electrical work, not of ionic currents. Integr. Physiol. Behav. Sci. 30, 283–306. doi:10.1007/BF02691602

PubMed Abstract | CrossRef Full Text | Google Scholar

Vrbjar, N., Dzurba, A., and Ziegelhoffer, A. (1994). Enzyme kinetics and the activation energy of Na/K-ATPase in ischaemic hearts: influence of the duration of ischaemia. Gen. Physiol. Biophys. 13, 405–411.

Google Scholar

Wang, L., Zhou, P., Craig, R. W., and Lu, L. (1999). Protection from cell death by mcl-1 is mediated by membrane hyperpolarization induced by K(+) channel activation. J. Membr. Biol. 172, 113–120. doi:10.1007/s002329900589

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, M., and Kondo, S. (2015). Is pigment patterning in fish skin determined by the turing mechanism? Trends Genet. 31, 88–96. doi:10.1016/j.tig.2014.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Werner, S., Stückemann, T., Beirá Amigo, M., Rink, J. C., Jucher, F., and Friedrich, B. M. (2015). Scaling and regeneration of self-organized patterns. Phys. Rev. Lett. 114, 138101. doi:10.1103/PhysRevLett.114.138101

PubMed Abstract | CrossRef Full Text | Google Scholar

Wright, S. (2004). Generation of resting membrane potential. Adv. Physiol. Educ. 28, 139–142. doi:10.1152/advan.00029.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamashita, M. (2013). Electric axon guidance in embryonic retina: galvanotropism revisited. Biochem. Biophys. Res. Commun. 431, 280–283. doi:10.1016/j.bbrc.2012.12.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M. (2009). Electrical fields in wound healing: an overriding signal that directs cell migration. Semin. Cell Dev. Biol. 20, 674–682. doi:10.1016/j.semcdb.2008.12.009

CrossRef Full Text | Google Scholar

Zhou, Y., Wei, X., Lu, M., Deng, B., Wang, J., and Che, Y. (2012). “A model of the endogenous electrical field effect: can ephaptic transmission cause neural synchrony independently?,” in Proceedings of the 31st Chinese Control Conference, Hefei.

Google Scholar