Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine

Simulations 1, 2, and 3 were used to validate the core BETSE model by determining its ability to predict resting Vmem and expected Vmem 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.

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

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 Vmem 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 Vmem and intracellular ion concentrations calculated by BETSE showed <10% difference between experimentally determined values measured from Xenopus oocytes (Table 2).

Simulation 2 explored resting Vmem 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.

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 (Vmem = 0 in all cells).

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

The simulation shows that after 20 simulated minutes, the BETSE-calculated Vmem 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 Vmem 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 Vmem. The Nernst equation [equation (21)] was used with BETSE-simulated intra- and extracellular dye concentrations as an alternative Vmem estimate (“Vmem Dye,” Table 3); results are virtually identical between the direct-BETSE and dye-estimated Vmem 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 Vmem 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 Vmem hyperpolarization (Lodish et al., 2000; Wright, 2004).

Simulation 3 explored factors influencing resting Vmem in isolated cells, and also demonstrated the ability for cell Vmem 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 Vmem (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 Vmem 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 Vmem, as listed for the 20 min time point in Table 3.

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

“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 Vmem is calculated assuming the cell is a simple capacitor via the charge inside the cell and the relation Vmem=1CmQcell. 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 Vmem 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 Vmem. 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 Ca2+ 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 Vmem for both simulation modes converged at precisely the Goldman Vmem 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 Vmem depolarization (Delamere and Duncan, 1977).

As can be seen in Figure 6, cell Vmem naturally returns to its resting value of −63.7 mV after each perturbation is complete. While overall, Vmem 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).

Simulation set 4 validated the expected function of voltage-gated Na+ and K+ channels, highlighted the ability for resting Vmem to control cell excitability, and examined the possibility of low voltage-gated K + expression in relation to voltage-gated Na + expression to effect resting Vmem. 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 Vmem, 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 Vmem 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.

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 Vmem 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 Vmem 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 Vmem (Figure 7C), which may assist in explaining the sustained depolarization of some cancer cells (Djamgoz, 2014).

Simulation 5 investigated the physiological impacts of a heterogeneous Vmem 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−5, and GJ connectivity with a value of βGJo = 1 × 10−7was 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 Vmem featuring a depolarized spot of cells in the lower left side (Vmem ~ −20 mV) and a hyperpolarized spot of cells in the upper right side of a circular cell cluster (Vmem ~ −60 mV). Each region was surrounded by cells with a mean Vmem of approximately −45 mV (Figure 8A).

The presence of regional Vmem differences was found to have various influences on the cluster as a whole. Heterogeneous Vmem induced differences in cytosolic Ca2+ levels (Figure 8B) in a manner inversely proportional to cell Vmem, with the most hyperpolarized cells having cytosolic Ca2+ 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 Vmem values (inverse to that of cationic Ca2+), 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 Vmem was also seen to produce significant osmotic pressure differences between cells of different resting potential, with more hyperpolarized Vmem leading to lower osmotic pressure than more depolarized Vmem (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 Vmem, closing to minimum conductivity values and isolating the two regions of differential Vmem from the remainder of the cluster (Figure 8D). The Vmem pattern in this example generated electric fields of up to ~6.5 × 105 V/m acting (over short spatial distances of 26 nm) between interfacing cell membranes of GJ networked cells (Figure 8G).

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

We conclude that stable patterns of resting Vmem 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.

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

Cells in a first simulation (low GJ connectivity) had an intercellular (GJ mediated) free-diffusion scaling of βGJo = 1.0 × 10−7. For a cell with 1.0 × 105 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 βGJo = 1.0 × 10−6. For a cell with 1.0 × 105 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 Vmem 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 Vmem, which remains stable as a new resting Vmem state divergent from that of surrounding cells (Figure 9D). However, the cluster with high GJ connectivity shows only 10 mV depolarization in Vmem 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 Vmem 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 Vmem states from a homogeneous collective.