Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine
Bioelectric cell properties have been revealed as powerful targets for modulating stem cell function, regenerative response, developmental patterning, and tumor reprograming. Spatio-temporal distributions of endogenous resting potential, ion flows, and electric fields are influenced not only by the genome and external signals but also by their own intrinsic dynamics. Ion channels and electrical synapses (gap junctions) both determine, and are themselves gated by, cellular resting potential. Thus, the origin and progression of bioelectric patterns in multicellular tissues is complex, which hampers the rational control of voltage distributions for biomedical interventions. To improve understanding of these dynamics and facilitate the development of bioelectric pattern control strategies, we developed the BioElectric Tissue Simulation Engine (BETSE), a finite volume method multiphysics simulator, which predicts bioelectric patterns and their spatio-temporal dynamics by modeling ion channel and gap junction activity and tracking changes to the fundamental property of ion concentration. We validate performance of the simulator by matching experimentally obtained data on membrane permeability, ion concentration and resting potential to simulated values, and by demonstrating the expected outcomes for a range of well-known cases, such as predicting the correct transmembrane voltage changes for perturbation of single cell membrane states and environmental ion concentrations, in addition to the development of realistic transepithelial potentials and bioelectric wounding signals. In silico experiments reveal factors influencing transmembrane potential are significantly different in gap junction-networked cell clusters with tight junctions, and identify non-linear feedback mechanisms capable of generating strong, emergent, cluster-wide resting potential gradients.
The BETSE platform will enable a deep understanding of local and long-range bioelectrical dynamics in tissues, and assist the development of specific interventions to achieve greater control of pattern during morphogenesis and remodeling.
Keywords: bioelectric simulation, pattern formation, resting potential, transmembrane voltage
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 (Vmem) 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 Vmem 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).
On the scale of single cells, the Vmem 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 Vmem (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 Vmem 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 Vmem 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 Vmem 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.
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.
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−, Ca2+, HCO3−, 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 Yn− or Yn+, 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 Vmem 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 Ca2+ 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.
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
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).
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 (ci) and intracellular voltage (Vcell) are defined. Each cell also has a unique volume (volcell) and perimeter representing the cell membrane. This allows unique membrane properties, such as Vmem, 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).
gradient (∇sj), which calculates the degree of change in the spatial property over space at grid point j
divergence (∇⋅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 (∇2sj = ∇ ⋅ ∇sj), 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 ∇2Sj = cj then Sj = ∇−2cj.
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.
Ion transport in bioelectrical systems is influenced by gradients of both concentration (∇ci) 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 ci of an ion i with charge zi and diffusion constant Di:
here, u→ is a fluid flow [e.g., electroosmotic flow field (Andreev, 2013)], q is the electron charge constant, kb is the Boltzmann constant, and T is the temperature (see Table 1).