Example 5: Mixed Topologies (Scalar Hub and Vectorized Populations)

This example demonstrates the flexibility of MTKNeuralToolkit's symbolic interfaces. You can seamlessly mix Scalar and Vectorized topologies in the same network. Because vectorized compartments expose their states as symbolic arrays, you can directly index into them (e.g., pop.interfaces.V[1]) to wire up standard scalar synapses to specific elements of a population.

We create a "hub" neuron (Scalar) that interacts with an E/I population pair (Vectorized). The hub excites the first neuron of the E-pop, which inhibits the hub. The E-pop also projects densely to the I-pop using a weight matrix.

using MTKNeuralToolkit
using MTKNeuralToolkit.HodgkinHuxley: SodiumChannel, PotassiumChannel, LeakChannel
using ModelingToolkit: mtkcompile, @named
using ModelingToolkitStandardLibrary.Blocks: Sine
using OrdinaryDiffEq
using Plots

1. Build the Scalar "Hub" Neuron

top_scalar = Scalar()

function build_hh_neuron(name::Symbol)
    @named cap  = Capacitor(topology=top_scalar, C=1.0)
    @named na   = SodiumChannel(topology=top_scalar)
    @named k    = PotassiumChannel(topology=top_scalar)
    @named leak = LeakChannel(topology=top_scalar)
    return build_compartment(cap, [na, k, leak]; name=name, V_init=-65.0, topology=top_scalar)
end

hub_neuron = build_hh_neuron(:hub_neuron)

Compartment{NoMorphology, NoGeometry, Float64}(Model hub_neuron:
Subsystems (8): see hierarchy(hub_neuron)
  cap
  injector
  syn_injector
  p
  ⋮
Equations (52):
  36 standard: see equations(hub_neuron)
  16 connecting: see equations(expand_connections(hub_neuron))
Unknowns (51): see unknowns(hub_neuron)
  cap₊v(t)
  cap₊i(t)
  cap₊p₊v(t)
  cap₊p₊i(t)
  ⋮
Parameters (7): see parameters(hub_neuron)
  cap₊C
  na₊g
  na₊E_rev
  k₊g
  ⋮, (V = hub_neuron₊cap₊v(t), p_pin = Model hub_neuron₊p:
Equations (2):
  2 connecting: see equations(expand_connections(hub_neuron₊p))
Unknowns (2): see unknowns(hub_neuron₊p)
  v(t)
  i(t), n_pin = Model hub_neuron₊n:
Equations (2):
  2 connecting: see equations(expand_connections(hub_neuron₊n))
Unknowns (2): see unknowns(hub_neuron₊n)
  v(t)
  i(t), I_ext = hub_neuron₊injector₊I₊u(t), I_syn = hub_neuron₊syn_injector₊I₊u(t), cap_name = :cap), -65.0, Scalar(), NoGeometry(), NoMorphology())

2. Build the Vectorized Populations (E and I)

N_E = 5
N_I = 3
top_E = Vectorized(N_E)
top_I = Vectorized(N_I)

function build_vector_pop(name::Symbol, top)
    @named cap  = Capacitor(topology=top, C=1.0)
    @named na   = SodiumChannel(topology=top)
    @named k    = PotassiumChannel(topology=top)
    @named leak = LeakChannel(topology=top)
    return build_compartment(cap, [na, k, leak]; name=name, V_init=-65.0, topology=top)
end

pop_E = build_vector_pop(:pop_E, top_E)
pop_I = build_vector_pop(:pop_I, top_I)

Compartment{NoMorphology, NoGeometry, Float64}(Model pop_I:
Subsystems (8): see hierarchy(pop_I)
  cap
  injector
  syn_injector
  p
  ⋮
Equations (64):
  48 standard: see equations(pop_I)
  16 connecting: see equations(expand_connections(pop_I))
Unknowns (51): see unknowns(pop_I)
  cap₊v(t)
  cap₊i(t)
  cap₊p₊v(t)
  cap₊p₊i(t)
  ⋮
Parameters (7): see parameters(pop_I)
  cap₊C
  na₊g
  na₊E_rev
  k₊g
  ⋮, (V = (pop_I₊cap₊v(t))[1:3], p_pin = Model pop_I₊p:
Equations (2):
  2 connecting: see equations(expand_connections(pop_I₊p))
Unknowns (2): see unknowns(pop_I₊p)
  v(t)
  i(t), n_pin = Model pop_I₊n:
Equations (2):
  2 connecting: see equations(expand_connections(pop_I₊n))
Unknowns (2): see unknowns(pop_I₊n)
  v(t)
  i(t), I_ext = (pop_I₊injector₊I₊u(t))[1:3], I_syn = (pop_I₊syn_injector₊I₊u(t))[1:3], cap_name = :cap), -65.0, Vectorized(3), NoGeometry(), NoMorphology())

3. Wire Synapses

A. Mixed Topology (Scalar <-> Vectorized)

The scalar hub excites the FIRST neuron in popE, and popE[1] inhibits the hub.

@named syn_hub_to_E = ExpSynapse(g_max=3.0, τ=5.0, E_rev=0.0, V_th=-20.0, slope=2.0)
@named syn_E_to_hub = ExpSynapse(g_max=5.0, τ=5.0, E_rev=-80.0, V_th=-20.0, slope=2.0)

Model syn_E_to_hub:
Equations (2):
  2 standard: see equations(syn_E_to_hub)
Unknowns (4): see unknowns(syn_E_to_hub)
  s(t)
  I_syn(t)
  V_pre(t)
  V_post(t)
Parameters (5): see parameters(syn_E_to_hub)
  g_max
  τ
  E_rev
  V_th
  ⋮

B. Vectorized <-> Vectorized

popE projects densely to popI using a weight matrix

W_EI = 0.5 .* rand(N_I, N_E)
syn_EI = build_synapse_block(pop_E, pop_I, W_EI; name=:syn_EI, E_rev=0.0)

synapse_specs = [
    #Hub -> pop_E[1]
    SynapseSpec(hub_neuron.interfaces.V, pop_E.interfaces.V[1], pop_E.interfaces.I_syn[1], syn_hub_to_E),
    # pop_E[1] -> Hub
    SynapseSpec(pop_E.interfaces.V[1], hub_neuron.interfaces.V, hub_neuron.interfaces.I_syn, syn_E_to_hub),
    #pop_E -> pop_I (Vectorized block)
    syn_EI
]

3-element Vector{SynapseSpec}:
 SynapseSpec(hub_neuron₊cap₊v(t), (pop_E₊cap₊v(t))[1], (pop_E₊syn_injector₊I₊u(t))[1], Model syn_hub_to_E:
Equations (2):
  2 standard: see equations(syn_hub_to_E)
Unknowns (4): see unknowns(syn_hub_to_E)
  s(t)
  I_syn(t)
  V_pre(t)
  V_post(t)
Parameters (5): see parameters(syn_hub_to_E)
  g_max
  τ
  E_rev
  V_th
  ⋮, nothing)
 SynapseSpec((pop_E₊cap₊v(t))[1], hub_neuron₊cap₊v(t), hub_neuron₊syn_injector₊I₊u(t), Model syn_E_to_hub:
Equations (2):
  2 standard: see equations(syn_E_to_hub)
Unknowns (4): see unknowns(syn_E_to_hub)
  s(t)
  I_syn(t)
  V_pre(t)
  V_post(t)
Parameters (5): see parameters(syn_E_to_hub)
  g_max
  τ
  E_rev
  V_th
  ⋮, nothing)
 SynapseSpec((pop_E₊cap₊v(t))[1:5], (pop_I₊cap₊v(t))[1:3], (pop_I₊syn_injector₊I₊u(t))[1:3], Model syn_EI:
Equations (2):
  2 standard: see equations(syn_EI)
Unknowns (4): see unknowns(syn_EI)
  s(t)
  I_syn(t)
  V_pre(t)
  V_post(t)
Parameters (6): see parameters(syn_EI)
  g_max
  τ
  E_rev
  V_th
  ⋮, Compartment{NoMorphology, NoGeometry, Float64}(Model pop_I:
Subsystems (8): see hierarchy(pop_I)
  cap
  injector
  syn_injector
  p
  ⋮
Equations (64):
  48 standard: see equations(pop_I)
  16 connecting: see equations(expand_connections(pop_I))
Unknowns (51): see unknowns(pop_I)
  cap₊v(t)
  cap₊i(t)
  cap₊p₊v(t)
  cap₊p₊i(t)
  ⋮
Parameters (7): see parameters(pop_I)
  cap₊C
  na₊g
  na₊E_rev
  k₊g
  ⋮, (V = (pop_I₊cap₊v(t))[1:3], p_pin = Model pop_I₊p:
Equations (2):
  2 connecting: see equations(expand_connections(pop_I₊p))
Unknowns (2): see unknowns(pop_I₊p)
  v(t)
  i(t), n_pin = Model pop_I₊n:
Equations (2):
  2 connecting: see equations(expand_connections(pop_I₊n))
Unknowns (2): see unknowns(pop_I₊n)
  v(t)
  i(t), I_ext = (pop_I₊injector₊I₊u(t))[1:3], I_syn = (pop_I₊syn_injector₊I₊u(t))[1:3], cap_name = :cap), -65.0, Vectorized(3), NoGeometry(), NoMorphology()))

4. Driving Stimuli & Network Assembly

Drive the hub neuron with a steady current to keep it spiking

@named current_driver = Sine(amplitude=8.0, frequency=0.05, offset=8.0)
drivers = [(1, current_driver)]

1-element Vector{Tuple{Int64, ModelingToolkitBase.System}}:
 (1, Model current_driver:
Subsystems (1): see hierarchy(current_driver)
  output
Equations (1):
  1 standard: see equations(current_driver)
Unknowns (1): see unknowns(current_driver)
  output₊u(t): Inner variable in RealOutput output
Parameters (5): see parameters(current_driver)
  offset
  start_time
  amplitude
  frequency
  ⋮)

Automatic Grounding

The network builder handles all grounding implicitly. Undriven elements in popE and popI will be grounded automatically.

net = build_acausal_network([hub_neuron, pop_E, pop_I];
                            synapse_specs=synapse_specs,
                            drivers=drivers,
                            name=:mixed_net)

println("Compiling mixed-topology network...")
sys = mtkcompile(net.sys)

Model mixed_net:
Equations (43):
  43 standard: see equations(mixed_net)
Unknowns (43): see unknowns(mixed_net)
  syn_hub_to_E₊s(t)
  syn_E_to_hub₊s(t)
  (syn_EI₊s(t))[5]
  (syn_EI₊s(t))[4]
  ⋮
Parameters (42): see parameters(mixed_net)
  hub_neuron₊cap₊C
  hub_neuron₊na₊g
  hub_neuron₊na₊E_rev
  hub_neuron₊k₊g
  ⋮
Observed (560): see observed(mixed_net)

Mixed systems generate block-diagonal Jacobians, so sparse solvers excel here.

prob = ODEProblem(sys, [], (0.0, 200.0), jac=true, sparse=true)

println("Solving...")
sol = solve(prob, Rosenbrock23())

retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 1690-element Vector{Float64}:
   0.0
   0.004151497658199239
   0.03346592387110653
   0.06681070233882805
   0.13339904295393767
   0.20698292714549937
   0.32143095284398615
   0.4475346744252218
   0.6156580192670218
   0.7978569217350175
   ⋮
 196.04154500837254
 196.68326857164058
 197.25722365095095
 197.8311787302613
 198.31690475171192
 198.80263077316252
 199.28835679461312
 199.71355563270956
 200.0
u: 1690-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.052, 0.052, 0.052  …  0.317, -65.0, -65.0, -65.0, -65.0, -65.0, 0.052, 0.596, 0.317, -65.0]
 [7.079314749575899e-13, 7.020810092988204e-13, 7.020810092987783e-13, 7.020810092987783e-13, 7.020810092987783e-13, 7.020810092987783e-13, 7.020810092988204e-13, 0.052016202304240935, 0.052016202304240935, 0.052016202304240935  …  0.3170005140127964, -65.00010510490591, -65.00010510490591, -65.00010510490591, -65.00010510490591, -65.00010510490563, 0.05201799239933115, 0.5959997767806697, 0.31700070751311876, -64.96691768807783]
 [6.032653711705216e-12, 5.642129272287841e-12, 5.642129272266868e-12, 5.642129272266868e-12, 5.642129272266868e-12, 5.642129272266868e-12, 5.642129272287841e-12, 0.052122661978541596, 0.05212266197854247, 0.05212266197854257  …  0.3170041020566747, -65.000722015743, -65.000722015743, -65.000722015743, -65.000722015743, -65.00072201572375, 0.052233439405473904, 0.5959821805016563, 0.3170166308056937, -64.7345045564028]
 [1.2849045990017192e-11, 1.1224759728148542e-11, 1.1224759727994339e-11, 1.1224759727994339e-11, 1.1224759727994339e-11, 1.1224759727994339e-11, 1.1224759728148542e-11, 0.05222833946975759, 0.05222833946976369, 0.05222833946976433  …  0.31700810972153937, -65.00118131994385, -65.00118131994385, -65.00118131994385, -65.00118131994385, -65.00118131986449, 0.0526513704369107, 0.5959280499523424, 0.3170579444512221, -64.47241526075793]
 [2.9267006779477026e-11, 2.2260114732444684e-11, 2.2260114731162657e-11, 2.2260114731162657e-11, 2.2260114731162657e-11, 2.2260114731162657e-11, 2.2260114732444684e-11, 0.052399375523972815, 0.05239937552401692, 0.05239937552402142  …  0.31701596603952875, -65.00147553586197, -65.00147553586197, -65.00147553586197, -65.00147553586197, -65.00147553551984, 0.05394459224559858, 0.5957115524098755, 0.31721372121543695, -63.95474763876269]
 [5.2855388230239835e-11, 3.428550199177494e-11, 3.428550198675777e-11, 3.428550198675777e-11, 3.428550198675777e-11, 3.428550198675777e-11, 3.428550199177494e-11, 0.052540127169883676, 0.05254012717003321, 0.05254012717004848  …  0.3170245503200039, -65.00106115462009, -65.00106115462009, -65.00106115462009, -65.00106115462009, -65.00106115371952, 0.05595605927155707, 0.5953035565919363, 0.3174989196843249, -63.389006781198]
 [1.0498813363859912e-10, 5.2650868497853855e-11, 5.265086847696368e-11, 5.265086847696368e-11, 5.265086847696368e-11, 5.265086847696368e-11, 5.2650868497853855e-11, 0.05269110421325363, 0.052691104213745794, 0.05269110421379603  …  0.3170379976101585, -64.99941883257983, -64.99941883257983, -64.99941883257983, -64.99941883257983, -64.9994188300485, 0.060006658773802554, 0.5943134190494644, 0.318177337634808, -62.515038386845966]
 [1.9532434788869737e-10, 7.242850165068308e-11, 7.242850158679555e-11, 7.242850158679555e-11, 7.242850158679555e-11, 7.242850158679555e-11, 7.242850165068308e-11, 0.05279671798126054, 0.052796717982425474, 0.05279671798254438  …  0.3170532451708352, -64.99679081354843, -64.99679081354843, -64.99679081354843, -64.99679081354843, -64.99679080765422, 0.0655110031555538, 0.5927075647024893, 0.31925826386324496, -61.548736464723376]
 [4.072633703272217e-10, 9.807992767347972e-11, 9.80799274711402e-11, 9.80799274711402e-11, 9.80799274711402e-11, 9.80799274711402e-11, 9.807992767347972e-11, 0.052884026707233, 0.05288402670981668, 0.05288402671008038  …  0.31707459639492674, -64.99272376175328, -64.99272376175328, -64.99272376175328, -64.99272376175328, -64.99272374706273, 0.07429265816783061, 0.5896821204652163, 0.3212563829852874, -60.22784662925493]
 [8.731969932385329e-10, 1.249881031310638e-10, 1.249881025704803e-10, 1.249881025704803e-10, 1.249881025704803e-10, 1.249881025704803e-10, 1.249881031310638e-10, 0.0529435675033784, 0.05294356750817578, 0.052943567508665475  …  0.31709918431099626, -64.98816708650033, -64.98816708650033, -64.98816708650033, -64.98816708650033, -64.98816705184163, 0.08576757908003876, 0.5851469169192259, 0.3241849676792076, -58.70098535555623]
 ⋮
 [0.10745351534949783, 0.12427395733283345, 8.460663727585066e-10, 8.460663727585068e-10, 8.460663727585066e-10, 8.460663727585068e-10, 0.12427395733283345, 0.02417228915245752, 0.024783592656554353, 0.04952039208796029  …  0.4544721427747331, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -60.806327913905854, 0.01947546388570156, 0.7036501791051311, 0.2573795107103382, -72.99933689697939]
 [0.09450256845477946, 0.10929570988781732, 8.460663727585105e-10, 8.460663727585107e-10, 8.460663727585105e-10, 8.460663727585107e-10, 0.10929570988781732, 0.026457207386664558, 0.027114441730602786, 0.050844975049045886  …  0.44585256919979294, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -61.21729846133239, 0.02140067157697537, 0.712479683118222, 0.25239115150054037, -72.1754876522877]
 [0.08424878281443271, 0.09743682873127357, 8.460663727585135e-10, 8.460663727585137e-10, 8.460663727585135e-10, 8.460663727585137e-10, 0.09743682873127357, 0.02864158443494514, 0.02933002757197054, 0.05179173494243401  …  0.4384021267572774, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -61.56391082908989, 0.024001031541854605, 0.7182371850074881, 0.24948792977831327, -71.17668672780488]
 [0.07510756080007053, 0.08686466834006328, 8.460663727585161e-10, 8.460663727585164e-10, 8.460663727585161e-10, 8.460663727585164e-10, 0.08686466834006328, 0.030917611549264958, 0.031626850836400704, 0.05252782146462519  …  0.4312012441583582, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -61.88959607245457, 0.027741602135044178, 0.7216091880552067, 0.24829699099452984, -69.90885252534193]
 [0.06815186981683302, 0.07882015512968633, 8.460663727585182e-10, 8.460663727585184e-10, 8.460663727585182e-10, 8.460663727585184e-10, 0.07882015512968633, 0.032897564606922995, 0.033615748248227884, 0.052999899007583305  …  0.4253047681233883, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -62.15066744468296, 0.03216593115800597, 0.7222714094582693, 0.24883654635835473, -68.60704097294088]
 [0.061840343513137205, 0.07152064208135656, 8.460663727585201e-10, 8.460663727585203e-10, 8.460663727585201e-10, 8.460663727585203e-10, 0.07152064208135656, 0.03489616106190376, 0.035615469325733645, 0.05334754920264861  …  0.4195910724754108, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -62.39749477472165, 0.03813390382301771, 0.720553283731378, 0.2510156486011712, -67.08491732772423]
 [0.0561133259857955, 0.06489713493737345, 8.460663727585217e-10, 8.460663727585219e-10, 8.460663727585217e-10, 8.460663727585219e-10, 0.06489713493737345, 0.036895490185420786, 0.037608202473067, 0.05358531086340534  …  0.4140617474579532, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -62.63079547863046, 0.046265780755328806, 0.7159894640491349, 0.25508780848736945, -65.3185379760539]
 [0.0515374372670154, 0.059604950659941904, 8.460663727585231e-10, 8.460663727585233e-10, 8.460663727585231e-10, 8.460663727585233e-10, 0.059604950659941904, 0.038632220816506385, 0.03933297748247413, 0.05371493225864994  …  0.4093731753743455, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -62.82461208250617, 0.05598678227128186, 0.709135619453282, 0.2604694709500175, -63.52482469769627]
 [0.048667524754822666, 0.056285790762638455, 8.460663727585239e-10, 8.460663727585241e-10, 8.460663727585239e-10, 8.460663727585241e-10, 0.056285790762638455, 0.03978870629897472, 0.04047825142650692, 0.05376624207472528  …  0.4062943656647368, -64.9997224337346, -64.9997224337346, -64.9997224337346, -64.9997224337346, -62.94994868347158, 0.06457115490155503, 0.7026568494961192, 0.26522183280110667, -62.135103811692694]

5. Plot the Results

p1 = plot(sol, idxs=[sys.hub_neuron.cap.v],
          title="Scalar Hub Neuron", ylabel="V (mV)", legend=false)

p2 = plot(sol, idxs=[sys.pop_E.cap.v...],
          title="Excitatory Pop (Pop[1] synapsed)", ylabel="V (mV)", legend=false)

p3 = plot(sol, idxs=[sys.pop_I.cap.v...],
          title="Inhibitory Pop (Driven by E)", ylabel="V (mV)", xlabel="Time (ms)", legend=false)

plot(p1, p2, p3, layout=(3,1), size=(800, 600))

This page was generated using Literate.jl.