Compartments

Single Compartments

A CompartmentSystem is sufficient to model a simple neuron–either as an isopotential spherical neuron, or as a "point" neuron with no morphological geometry. At minimum, a CompartmentSystem will contain information about its intrinsic membrane currents and ionic gradients (i.e. equilibrium potentials). We define these by directly providing vectors of ConductanceSystem and EquilibriumPotential to the CompartmentSystem:

using Conductor, Unitful, ModelingToolkit
import Unitful: mV, mS, cm

Vₘ = ParentScope(MembranePotential())

nav_kinetics = [
    Gate(AlphaBeta,
         ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         4.0*exp(-(Vₘ + 65.0)/18.0), p = 3, name = :m)
    Gate(AlphaBeta,
         0.07*exp(-(Vₘ+65.0)/20.0),
         1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)), name = :h)]

kdr_kinetics = [
    Gate(AlphaBeta,
         ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         0.125 * exp(-(Vₘ + 65.0)/80.0),
         p = 4, name = :n)]

# Note: `IonChannel` is a convenience constructor that returns a `ConductanceSystem`
@named NaV = IonChannel(Sodium, nav_kinetics, max_g = 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, max_g = 36mS/cm^2)
@named leak = IonChannel(Leak, max_g = 0.3mS/cm^2)

channels = [NaV, Kdr, leak];
reversals = Equilibria([Sodium => 50.0mV, Potassium => -77.0mV, Leak => -54.4mV])
dynamics = HodgkinHuxley(channels, reversals)
@named neuron = CompartmentSystem(Vₘ, dynamics)

\[ \begin{align} \frac{\mathrm{d} V_m\left( t \right)}{\mathrm{d}t} =& \frac{ - Kdr_{+}I\left( t \right) - NaV_{+}I\left( t \right) - leak_{+}I\left( t \right)}{a_m c_m} \\ \frac{\mathrm{d} NaV_{+}m\left( t \right)}{\mathrm{d}t} =& \left( 1 - NaV_{+}m\left( t \right) \right) \mathrm{ifelse}\left( 1, \frac{0.1 \left( 40 + V_m\left( t \right) \right)}{1 - e^{0.1 \left( -40 - V_m\left( t \right) \right)}}; \left( V_m\left( t \right) = -40 \right) \right) - 4 NaV_{+}m\left( t \right) e^{0.055556 \left( -65 - V_m\left( t \right) \right)} \\ \frac{\mathrm{d} NaV_{+}h\left( t \right)}{\mathrm{d}t} =& \frac{ - NaV_{+}h\left( t \right)}{1 + e^{0.1 \left( -35 - V_m\left( t \right) \right)}} + 0.07 \left( 1 - NaV_{+}h\left( t \right) \right) e^{0.05 \left( -65 - V_m\left( t \right) \right)} \\ NaV_{+}g\left( t \right) =& \left( NaV_{+}m\left( t \right) \right)^{3} NaV_{+}gbar NaV_{+}h\left( t \right) \\ NaV_{+}I\left( t \right) =& a_m \left( - ENa + V_m\left( t \right) \right) NaV_{+}g\left( t \right) \\ \frac{\mathrm{d} Kdr_{+}n\left( t \right)}{\mathrm{d}t} =& \left( 1 - Kdr_{+}n\left( t \right) \right) \mathrm{ifelse}\left( 0.1, \frac{0.01 \left( 55 + V_m\left( t \right) \right)}{1 - e^{0.1 \left( -55 - V_m\left( t \right) \right)}}; \left( V_m\left( t \right) = -55 \right) \right) - 0.125 Kdr_{+}n\left( t \right) e^{0.0125 \left( -65 - V_m\left( t \right) \right)} \\ Kdr_{+}g\left( t \right) =& \left( Kdr_{+}n\left( t \right) \right)^{4} Kdr_{+}gbar \\ Kdr_{+}I\left( t \right) =& a_m \left( - EK + V_m\left( t \right) \right) Kdr_{+}g\left( t \right) \\ leak_{+}g\left( t \right) =& leak_{+}gbar \\ leak_{+}I\left( t \right) =& a_m \left( - El + V_m\left( t \right) \right) leak_{+}g\left( t \right) \end{align} \]

In the case above, the CompartmentSystem constructor assumes a dimensionless geometry = Point(). The maximum magnitudes of the ion channel conductances, $\overline{g}$, have units of SpecificConductance (mS/cm²) that scale with the surface area of the compartment. However, when the geometry is a Point(), the CompartmentSystem ignores surface area scaling (internally, the area is fixed to 1.0). We can instead provide a <:Geometry object to describe the shape and size of the compartment:

import Unitful: µm
soma_shape = Sphere(radius = 20µm)
geo_dynamics = HodgkinHuxley(channels, reversals)
@named neuron = CompartmentSystem(Vₘ, geo_dynamics; geometry = soma_shape)

If we try to simulate the neuron we've modeled so far, the result isn't very interesting:

import Unitful: ms
using OrdinaryDiffEq, Plots
sim = Simulation(neuron, 300ms)
solution = solve(sim, Rosenbrock23())
plot(solution; idxs=[Vₘ])

The neuron isn't spontaneously active. To make the neuron produce spikes, we can write an equation for an electrode current and provide it to CompartmentSystem:

import Unitful: µA, pA

# A 400 picoamp squarewave pulse when 100ms > t > 200ms
@named step_pulse = PulseTrain(amplitude = 400.0pA, duration = 100ms, delay = 100ms)

stim_dynamics = HodgkinHuxley(channels, reversals)

@named neuron_stim = CompartmentSystem(Vₘ, stim_dynamics;
                                       geometry = soma_shape,
                                       stimuli = [step_pulse])

Putting it all together, our neuron simulation now produces a train of action potentials:

sim = Simulation(neuron_stim, 300ms)
solution = solve(sim, Rosenbrock23())
plot(solution; idxs=[Vₘ])

Multiple Compartments

Neurons with multiple compartments are explicitly constructed by defining a MultiCompartmentTopology, a directed graph that is provided to MultiCompartmentSystem. Individual subcompartments are connected via add_junction!. By default, a junction between compartments is assumed to be symmetrical–the axial conductance from compartment A to B is the same value as from compartment B to A. However, if the conductance between compartments is asymmetric, two junctions must be defined–one for the conductance in each direction.

For example usage, see the pinskyrinzel.jl demo.

struct CompartmentSystem{T<:Conductor.AbstractDynamics} <: Conductor.AbstractCompartmentSystem

struct MultiCompartmentSystem <: Conductor.AbstractCompartmentSystem

MultiCompartment(
    topology::MultiCompartmentTopology;
    extensions,
    name,
    defaults
) -> MultiCompartmentSystem