EMI basics

See code in GitLab.
Author: Tobias Gerach tobias.gerach@kit.edu

EMI basics.

The EMI model (Extracellular–Membrane–Intracellular) resolves each cardiomyocyte as an independent intracellular subdomain coupled to a shared extracellular domain and, if present, to neighbouring myocytes through explicit membrane faces. This is fundamentally different from bidomain or monodomain, where all intracellular space is treated as a continuous medium.

This example introduces the building blocks of an EMI simulation on a minimal geometry: two 50 × 18 × 18 µm3 cardiomyocytes lying end to end in an extracellular domain, coupled by a gap junction at their shared face. It covers the EMI mesh anatomy, the EMI-specific carputils functions, gap-junction and per-cell heterogeneity, and how to visualise the extracellular potential \(\phi_e\) in meshalyzer.

EMI mesh setup

A standard openCARP mesh (.pts/.elem) describes a uniform volume. EMI additionally requires two tag-list files that declare which element tags form intracellular domains and which form the extracellular domain:

<mesh>.intra

One intracellular region per cardiomyocyte. Format: first line is the count, then one tag per line.

<mesh>.extra

Tags of the extracellular domain. Same format.

openCARP auto-detects the membrane interface as any face shared between an intracellular and an extracellular element. Gap junctions are faces shared between two intracellular elements.

The mesh is generated on the fly with mesh.Block and two mesh.BoxRegion tags, one per cell. BoxRegion has an asymmetric convention: lowerleft is the absolute lower corner, but the mesher adds the block centre to upperright, so the upper corner must be supplied relative to the centre. A small helper hides this and tags a cell by its absolute x-range:

centre = (90.5, 65, 13)   # µm
geom = mesh.Block(size=(0.181, 0.130, 0.026),   # mm: 181 × 130 × 26 µm
                  resolution=0.002, etype='tetra',
                  centre=tuple(c/1000 for c in centre))

def cell_box(xmin, xmax, tag):
    lo = (xmin, 56, 4)
    hi = (xmax - centre[0], 74 - centre[1], 22 - centre[2])
    return mesh.BoxRegion(lowerleft=tuple(v/1000 for v in lo),
                          upperright=tuple(v/1000 for v in hi),
                          tag=tag, bath=0)

geom.add_region(cell_box(40, 90,  tag=2))   # left cell
geom.add_region(cell_box(90, 140, tag=4))   # right cell

The two boxes abut at the midplane; the face they share becomes the gap junction.

After generation the script writes the two EMI tag files alongside the mesh:

# .extra: one region, tag 1 (the extracellular domain — default mesher tag)
with open(meshname + '.extra', 'w') as f:
    f.write('1\n1')

# .intra: two regions, tags 2 and 4 (the two cells)
with open(meshname + '.intra', 'w') as f:
    f.write('2\n2\n4')

Parallel execution

EMI partitions the mesh one tag at a time: every element sharing an intra- or extracellular tag is assigned to the same MPI rank, so a cell or the surrounding extracellular domain is never split across ranks. The number of MPI ranks is therefore bounded by the total number of unique tags — the intracellular tags listed in <mesh>.intra plus the extracellular tags in <mesh>.extra:

n_ranks <= (unique intracellular tags) + (unique extracellular tags)

Tags are distributed to ranks in ascending order by default; an optional <mesh>.part file (one tag / rank pair per line) overrides this mapping when a specific assignment is desired. Requesting more ranks than there are tags aborts the run with an explicit error rather than silently degrading the partition.

EMI sub-meshes in the output

During setup openCARP decouples the mesh into two sub-meshes written to the simulation output directory:

<basename>_e.*

Volumetric mesh of the decoupled domain — the extracellular region and both intracellular cells. phie.igb is stored here as nodal data: the extracellular potential \(\phi_e\) on extracellular nodes and the intracellular potential on the cell nodes. Being nodal, meshalyzer can plot time traces at individual nodes.

<basename>_m.*

Surface mesh of the membrane interface. vm.igb is stored here as cell data — one value per membrane face, not per node. Meshalyzer renders it as a colour map; it cannot plot per-element time traces.

Physics region

tools.gen_physics_opts registers the EMI physics type with openCARP:

tools.gen_physics_opts(EMI=True)

generates

-num_phys_regions 1
-phys_region[0].name "EMI domain"
-phys_region[0].ptype 6

ptype 6 is the EMI physics identifier. For bidomain/monodomain simulations ptype 0 (intracellular) and ptype 1 (extracellular) are used instead.

Conductivity regions

EMI uses model.ConductivityRegionEMI in place of the bidomain model.ConductivityRegion. The two factory methods extra_default and intra_default create catch-all regions that apply to all extracellular and intracellular elements respectively, without requiring an explicit tag list. openCARP determines which elements belong to each category from the .extra and .intra tag files.

model.ConductivityRegionEMI.extra_default(g_bath=2.0)
model.ConductivityRegionEMI.intra_default(g_bath=0.4)

These two calls, passed via opts_formatted(i), generate:

-num_gregions 2
-gregion[0].name "extra_regions_default"
-gregion[0].g_bath 2.0
-gregion[0].g_mult 1.0
-gregion[1].name "intra_regions_default"
-gregion[1].g_bath 0.4
-gregion[1].g_mult 1.0

g_bath is the isotropic conductivity in S/m. No longitudinal/transverse split is needed because EMI resolves individual cells and their geometry explicitly rather than relying on fibre-direction homogenisation.

Individual region conductivities can still be modified by tag ID to override the defaults for specific element groups. For example, to raise the intracellular conductivity of the right-hand cell (tag 4) relative to the default:

model.ConductivityRegionEMI([4], 'region 4', g_bath=1.5)

which generates:

-gregion[2].name "region 4"
-gregion[2].num_IDs 1
-gregion[2].ID[0] 4
-gregion[2].g_bath 1.5
-gregion[2].g_mult 1.0

The two default regions already cover both cells and the extracellular domain; an explicit entry like the one above is only needed to make an individual cell heterogeneous.

Stimulation

The cell pair is excited by a transmembrane current stimulus applied at the left end. Dirichlet grounds on the outer faces of the extracellular domain pin its potential here, but they are not required: with no ground or Dirichlet boundary the solver handles the floating extracellular potential through null-space removal. Stimulus modes and electrode geometry are covered in the 01_stimulation example.

Ionic model regions

EMI requires two ionic model regions, one for membrane faces (between an intracellular and an extracellular element) and one for gap-junction faces (between two intracellular elements). The emi=True keyword switches the carputils array from imp_region to imp_region_emi:

model.ionic.AlievPanfilovIonicModel(None, 'default_ionic', emi=True)
model.ionic.PlonseyIonicModel(None, 'default_gap_junction', Rm=0.0045, emi=True)

generates:

-num_imp_regions 2
-imp_region_emi[0].name "default_ionic"
-imp_region_emi[0].im AlievPanfilov
-imp_region_emi[1].name "default_gap_junction"
-imp_region_emi[1].im Plonsey
-imp_region_emi[1].im_param Rm=0.0045

openCARP assigns imp_region_emi[0] to membrane faces and imp_region_emi[1] to gap-junction faces based on the mesh topology. IDs=None (first positional argument) means the model applies to all faces of that type; no explicit face-tag list is needed.

The Plonsey model represents a gap junction as a purely resistive membrane with resistance \(R_m\) [kΩ cm²]. Here the two cells share a single internal face, so imp_region_emi[1] is active on that one gap junction; with the low default resistance the pair conducts almost like a continuous cell.

Heterogeneity can be introduced by adding extra imp_region_emi entries with explicit face IDs. In the EMI model, the gap-junction face between two cells with element tags a and b is identified by the compound face ID a:b. The junction between cells 2 and 4 is therefore 2:4; to raise its resistance and slow conduction across it:

model.ionic.PlonseyIonicModel(['2:4'], 'blocking_gj', Rm=0.5, emi=True)

This becomes imp_region_emi[2] — a third entry that overrides the default gap-junction region for that face group. Membrane faces are identified the same way, by the pair of the intracellular tag and the extracellular tag: cell 4's membrane is 4:1. Giving it its own imp_region_emi entry assigns that cell a distinct ionic model — the basis for representing different cell types (atrial, ventricular, nodal):

model.ionic.tenTusscherPanfilovIonicModel(['4:1'], 'tt_cell_4', emi=True)

Visualisation

With --visualize, meshalyzer opens phie.igb on the volumetric sub-mesh <basename>_e with a pre-configured view (clipping plane through the block midpoint, potential range ±0.0125 mV at time slice 10, isolines over the same range). This view is tuned for the default homogeneous mode; the other modes produce different ranges and need re-scaling. To inspect the transmembrane voltage open vm.igb on <basename>_m in a second window.

Usage

cd examples/07_EMI/00_basics
./run.py --visualize

--mode selects a heterogeneity demonstration on the right cell (tag 4): conductivity overrides its intracellular conductivity, block raises the gap-junction resistance to block conduction across the 2:4 interface, and ionic switches it to the tenTusscher-Panfilov model. The default homogeneous leaves both cells identical.

./run.py --mode block --visualize

Recent questions tagged EMI, carputils