# Two cells connected via a gap junction¶

In this example, we will set up two cells connected via a gap junction. The cells have different leak potentials. We will investigate how the equilibrium potentials of the two cells change because of the gap junction connection.

Note

Concepts covered in this example:

1. Creating a simulation recipe for two cells.

2. Placing probes.

3. Running the simulation and extracting the results.

4. Adding a gap junction connection.

## The full code¶

You can find the full code of the example at `python/examples/two_cell_gap_junctions.py`

Executing the script will run the simulation with default parameters.

## Walk-through¶

We set up a recipe for the simulation of two cells with some parameters.

```#!/usr/bin/env python3

import arbor
import argparse
import numpy as np

import pandas
import seaborn  # You may have to pip install these.
import matplotlib.pyplot as plt

class TwoCellsWithGapJunction(arbor.recipe):
def __init__(self, probes,
Vms, length, radius, cm, rL, g, gj_g,
cv_policy_max_extent):
"""
probes -- list of probes

Vms -- membrane leak potentials of the two cells
length -- length of cable in μm
radius -- radius of cable in μm
cm -- membrane capacitance in F/m^2
rL -- axial resistivity in Ω·cm
g -- membrane conductivity in S/cm^2
gj_g -- gap junction conductivity in μS

cv_policy_max_extent -- maximum extent of control volume in μm
"""

# The base C++ class constructor must be called first, to ensure that
# all memory in the C++ class is initialized correctly.
arbor.recipe.__init__(self)

self.the_probes = probes

self.Vms = Vms
self.length = length
self.cm = cm
self.rL = rL
self.g = g
self.gj_g = gj_g

self.cv_policy_max_extent = cv_policy_max_extent

self.the_props = arbor.neuron_cable_properties()

def num_cells(self):
return 2

def num_sources(self, gid):
assert gid in [0, 1]
return 0

def cell_kind(self, gid):
assert gid in [0, 1]
return arbor.cell_kind.cable

def probes(self, gid):
assert gid in [0, 1]
return self.the_probes

def global_properties(self, kind):
return self.the_props

def cell_description(self, gid):
```

Implementing the `cell_description` member function constructs the morphology and sets the properties of the cells as well as the gap junction mechanisms and the discretization policy.

```        return self.the_probes

def global_properties(self, kind):
return self.the_props

def cell_description(self, gid):
"""A high level description of the cell with global identifier gid.

For example the morphology, synapses and ion channels required
to build a multi-compartment neuron.
"""
assert gid in [0, 1]

tree = arbor.segment_tree()

tree.append(arbor.mnpos,
arbor.mpoint(0, 0, 0, self.radius),
arbor.mpoint(self.length, 0, 0, self.radius),
tag=1)

labels = arbor.label_dict({'cell' : '(tag 1)', 'gj_site': '(location 0 0.5)'})

decor = arbor.decor()
decor.set_property(Vm=self.Vms[gid])
decor.set_property(cm=self.cm)
decor.set_property(rL=self.rL)

# add a gap junction mechanism at the "gj_site" location and label that specific mechanism on that location "gj_label"
junction_mech = arbor.junction('gj', {"g" : self.gj_g})
decor.place('"gj_site"', junction_mech, 'gj_label')
decor.paint('"cell"', arbor.density(f'pas/e={self.Vms[gid]}', {'g': self.g}))

if self.cv_policy_max_extent is not None:
policy = arbor.cv_policy_max_extent(self.cv_policy_max_extent)
decor.discretization(policy)
else:
decor.discretization(arbor.cv_policy_single())

return arbor.cable_cell(tree, labels, decor)

def gap_junctions_on(self, gid):
```

The bidirection gap junction is created in the function `gap_junctions_on`.

```
# create a bidirectional gap junction from cell 0 at label "gj_label" to cell 1 at label "gj_label" and back.
return [arbor.gap_junction_connection((1 if gid == 0 else 0, 'gj_label'), 'gj_label', 1)]

```

We parse the command line arguments, instantiate the recipe, run the simulation, extract results and plot:

```    parser = argparse.ArgumentParser(description='Two cells connected via a gap junction')

'--Vms', help="membrane leak potentials in mV", type=float, default=[-100, -60], nargs=2)
'--length', help="cell length in μm", type=float, default=100)
'--radius', help="cell radius in μm", type=float, default=3)
'--cm', help="membrane capacitance in F/m^2", type=float, default=0.005)
'--rL', help="axial resistivity in Ω·cm", type=float, default=90)
'--g', help="membrane conductivity in S/cm^2", type=float, default=0.001)

'--gj_g', help="gap junction conductivity in μS", type=float, default=0.01)

help="maximum extent of control volume in μm", type=float)

# parse the command line arguments
args = parser.parse_args()

# set up membrane voltage probes at the position of the gap junction
probes = [arbor.cable_probe_membrane_voltage('"gj_site"')]
recipe = TwoCellsWithGapJunction(probes, **vars(args))

# create a default execution context and a default domain decomposition
context = arbor.context()
domains = arbor.partition_load_balance(recipe, context)

# configure the simulation and handles for the probes
sim = arbor.simulation(recipe, domains, context)

dt = 0.01
handles = []
for gid in [0, 1]:
handles += [sim.sample((gid, i), arbor.regular_schedule(dt)) for i in range(len(probes))]

# run the simulation for 5 ms
sim.run(tfinal=5, dt=dt)

# retrieve the sampled membrane voltages and convert to a pandas DataFrame
print("Plotting results ...")
df_list = []
for probe in range(len(handles)):
samples, meta = sim.samples(handles[probe])[0]
df_list.append(pandas.DataFrame({'t/ms': samples[:, 0], 'U/mV': samples[:, 1], 'Cell': f"{probe}"}))

df = pandas.concat(df_list,ignore_index=True)

fig, ax = plt.subplots()

# plot the membrane potentials of the two cells as function of time
seaborn.lineplot(ax=ax,data=df, x="t/ms", y="U/mV",hue="Cell",ci=None)

# area of cells
area = args.length*1e-6 * 2*np.pi*args.radius*1e-6

# total conductance and resistance
cell_g = args.g/1e-4 * area
cell_R = 1/cell_g

# gap junction conductance and resistance in base units
si_gj_g = args.gj_g*1e-6
si_gj_R = 1/si_gj_g

# indicate the expected equilibrium potentials
for (i, j) in [[0,1], [1,0]]:
weighted_potential = args.Vms[i] + ((args.Vms[j] - args.Vms[i])*(si_gj_R + cell_R))/(2*cell_R + si_gj_R)
ax.axhline(weighted_potential, linestyle='dashed', color='black', alpha=0.5)

# plot the initial/nominal resting potentials
for gid, Vm in enumerate(args.Vms):
ax.axhline(Vm, linestyle='dashed', color='black', alpha=0.5)

fig.savefig('two_cell_gap_junctions_result.svg')
```

The output plot below shows how the potential of the two cells approaches their equilibrium potentials. The expected values are denoted by dashed lines.

The equivalent circuit used to calculate the equilibrium potentials is depicted below.