bioelectrics

A differentiable engine for somatic bioelectricity, in PyTorch.

Documentation: https://tlcdv.github.io/bioelectrics/ — with four runnable tutorials that open in Colab with no installation — including a validation page whose every number is recomputed from this code, and a limitations page.

bioelectrics simulates the membrane potential (V_mem) of non-neural tissues at the level of ions, channels, pumps and gap junctions, and computes exact gradients of any outcome with respect to any parameter. It uses these gradients for the bioelectric inverse problem: which experimentally realizable intervention (a channel mRNA, a drug dose, a gap junction blocker, a light pattern), applied for a short time, moves a tissue into a target V_mem pattern that then persists on its own?

The library builds on the work of Michael Levin's lab, Alexis Pietak (BETSE), and Javier Cervera and Salvador Mafe. See Related tools for what it takes from them and what it adds.

(A) An ion-level cell (GHK fluxes, pumps, a K⁺ inward rectifier and a depolarizing cation channel) has two attractors. (B) A depolarized "tumor" domain is isolated by closed voltage-gated gap junctions. (C) An optimized bath drug, given for 0.1 s and then washed out, flips the tumor cells into the healthy attractor. (D) An optimized light pattern for an anion channelrhodopsin concentrates the light on the tumor. (E) The tissue after the light is off. (F) Optimization by backpropagation through the simulation.

Results of this demo (16×16 cells, 31 tumor cells):

The full protocol tables are in assets/ionic_protocols.md.

What you can do with it

Installation

git clone https://github.com/tlcdv/bioelectrics.git
cd bioelectrics
python -m venv .venv
.venv/bin/pip install -e ".[dev]"      # Windows: .venv\Scripts\pip
.venv/bin/python -m pytest             # about 3 minutes on a laptop CPU

Requires Python ≥ 3.10 and PyTorch ≥ 2.0. A GPU is not required.

Quickstart: design a drug protocol

import torch
from bioelectrics.ionic import IonicTissue, ChannelExpression, K_IR_CAM, CAT_DEP_CAM
from bioelectrics.interventions import ChannelBlocker, InterventionOptimizer

# Hardware: each cell expresses a K+ inward rectifier and a depolarizing cation channel.
tissue = IonicTissue(channels=[ChannelExpression(K_IR_CAM, 1e-16),
                               ChannelExpression(CAT_DEP_CAM, 1e-16)])

# Software: an 8x8 tissue at -78 mV with a 2x2 depolarized domain at -3 mV.
v0 = torch.full((1, 8, 8), -0.078, dtype=torch.float64)
v0[:, 3:5, 3:5] = -0.003
state = tissue.simulate(tissue.with_vmem(tissue.initial_state((8, 8)), v0), 0.5)

# Which dose of a bath-applied blocker of the depolarizing channel, given for 0.1 s and
# then washed out, puts every cell into the healthy attractor by t = 0.5 s?
opt = InterventionOptimizer(tissue, state, torch.full_like(v0, -0.078),
                            [ChannelBlocker("Cat_dep(CAM)", cost=1.0)], t_stim=0.1, t_relax=0.5)
print(opt.optimize(n_iters=30, tol_mV=1.0).protocol.to_markdown())

Examples

Concepts

The API uses the terms of DOMAIN_CONCEPTS.md, after Levin & Martyniuk (2018) and Levin (2019):

• Hardware and software. Channel and gap junction expression (ChannelExpression, IonicTissue parameters) is the hardware. The V_mem pattern (IonicState, BioelectricState) is the software, which stores pattern memory.
• Target morphology as an attractor. The optimizers are error-minimizing controllers that move a tissue from one basin of attraction into a target attractor.
• Triggers, not clamps. Every inverse design has an intervention phase [0, T_stim] and an intervention-free phase [T_stim, T_relax]. The error is measured at T_relax, so a solution must be a permanent attractor switch.
• Computational boundary of the self. A large transjunctional voltage closes gap junctions and isolates a depolarized domain; re-polarization re-opens them.

Physics and validation

Equations, parameters and their sources are in MATH_REFERENCE.md. All internal math is in SI units. The tests check the physics against analytic results and published results:

Related tools

Limitations and roadmap

Present limitations, stated so that results are not over-read:

• The extracellular space is a fixed bath (as BETSE without is_ecm). There are no extracellular electric fields or extracellular diffusion.
• The CAM-type channels used for bistability are phenomenological, not models of one specific protein. The BETSE Kir2.1 model alone does not give bistability with the default profile.
• The model has no cell volume or osmotic regulation, and the bath is fixed. A depolarized steady state is reached within seconds, but over hours (slowest time scale about 3 h) its cells gain Na⁺ (up to about 200 mM in the default cell). Steady states of depolarized cells are therefore not physiological at long times; see MATH_REFERENCE.md §7.
• bioelectrics.pharmacology contains 10 potency values read from 4 primary sources. Each carries its assay, species, regulatory status and citation. None of them was estimated. Compounds considered and left out are listed in pharmacology.EXCLUDED with the reason; dihydropyridines are excluded because their block is state-dependent, so one IC50 is wrong in a tissue whose cells sit at different voltages.
• The screen in examples/drug_repurposing.py finds that no drug in that set normalizes a depolarized domain. Every compound with a verified value acts on a polarizing element (the inward rectifier, the Na/K-ATPase), so blocking it depolarizes. This is a statement about which potency data are published, not about bioelectric medicine. The depolarizing channel in the model is also phenomenological rather than a named protein, so no real drug targets it by construction.
• Potency values are in-silico inputs, not predictions of a clinical effect. Off-target actions are recorded but not modelled.
• BPTT cost grows with simulated time; long rollouts on large ion-level grids are slow on a CPU (an 8×8 tissue: about 4 s per optimization iteration for a 0.5 s rollout).

Planned:

• BEN (Manicka & Levin 2019) in the zoo. The equations in its Fig. 8 have inconsistent units and signs, so an exact reproduction needs confirmation from the authors.
• Higher-order stiff solvers, GPU benchmarks, Voronoi cell meshes from segmented images.
• Gene regulatory network coupling (V_mem → expression → channels).
• Fitting model parameters to voltage-dye (DiBAC, CC2-DMPE) images.

Contributing and citation

See CONTRIBUTING.md. If you use the library, cite it with CITATION.cff, and cite the papers of the models you use.

References

• Levin, M. & Martyniuk, C. J. (2018). The bioelectric code: An ancient computational medium for dynamic control of growth and form. BioSystems 164, 76–93.
• Levin, M. (2019). The computational boundary of a "self": Developmental bioelectricity drives multicellularity and scale-free cognition. Frontiers in Psychology 10, 2688.
• Pietak, A. & Levin, M. (2016). Exploring instructive physiological signaling with the bioelectric tissue simulation engine. Frontiers in Bioengineering and Biotechnology 4, 55.
• Cervera, J., Alcaraz, A. & Mafe, S. (2016). Bioelectrical signals and ion channels in the modeling of multicellular patterns and cancer biophysics. Scientific Reports 6, 20403.
• Manicka, S. & Levin, M. (2019). Modeling somatic computation with non-neural bioelectric networks. Scientific Reports 9, 18612.
• Verwer, J. G., Spee, E. J., Blom, J. G. & Hundsdorfer, W. (1999). A second-order Rosenbrock method applied to photochemical dispersion problems. SIAM J. Sci. Comput. 20, 1456–1480.

License

MIT. Parts of src/bioelectrics/ionic/ are derived from BETSE (BSD-2-Clause); see THIRD_PARTY_NOTICES.md.