Exploring Instructive Physiological Signaling with the Bioelectric Tissue Simulation Engine

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 μAcm2 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.

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:

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 Vmem, which is used as a cross-check of BETSE-derived concentrations and Vmem calculations.

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” Vmem. Due to the action of active ion pumping in living cells, the steady-state Vmem represents a dynamic equilibrium with net ion flux but zero current, and can be estimated from the Goldman equation as:

In the Goldman equation (22), ions are separated into anions (c−) and cations (c+) with concentrations inside (cint) and outside (cext) of the cell membrane. The membrane has a specific permeability Pmem 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 Vmem 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 Vmem values obtained using BETSE’s Maxwell Capacitance Matrix voltage solving method defined in section 2.5.

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, Keqm, can be expressed both in terms of the reaction free energy under standard conditions, △Greacteqm, and in terms of the reaction’s product concentrations (index k) and those of its reactants (index j) where ak and aj represent coefficients of stoichiometry for the reaction (Beard and Qian, 2007; Pekar, 2015):

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

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 ak and aj represent coefficients of stoichiometry for the reaction:

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

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

when ΔGpump = 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 Vmem is:

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

Values for the Michaelis constants KNa = 5.0, KK = 0.2, and KATP = 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.

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

Gap junctions were modeled as (optionally) voltage-sensitive conduits influencing the intercellular diffusion coefficient for all ions uniformly via a diffusion–constant scaling factor, βGJo. 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 (βGJo=0). Medium-high GJ connectivity corresponded to βGJo = 1.0 × 10−6, an intercellular diffusion coefficient of approximately 1.0 × 10−15m2/s. Assuming 1.0 × 105 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 (VGJ) and time. Specific details regarding voltage gating of GJ are described in Supplementary Material.

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

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 rgj ~ 5 to 8 nm) and the narrow channels (decm ~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:

where F^e is a volume force generated by electrostatic forces resulting from a voltage gradient (electric field E→) between two cells or extracellular spaces:

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:

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

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

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.