This page was generated from docs/paper_tsodyksuzielmarkram2000.ipynb. Interactive online version: Colab badge Deepnote badge. Binder badge Download notebook.

Example: Tsodyks, Uziel, Markram (2000)

Here we see how to implement the network synchrony generation model from Tsodyks, Uziel, Markram (2000). This network exhibits network-wide synchronised firing (along with other activity), using only event-based synapses!

This example requires and demonstrates the use of two of EDEN’s extensions to NeuroML:

  • Parameter variability among point neurons and synapses;

  • Multiflux: Since the cell models the time constant rather than passive leak and capacitance of the membrane, stimuli come in terms of voltage! Such that \(\dot{V}_m = \frac{-V_m + I_{ext}}{\tau_m}\) . To be fair, in this case we could introduce physical leak and (arbitrary) capacitance to the model, and carefully adjust most parameters to the same current-based effect … but this starts being quite the burden, and we can’t always stay within NeuroML2 anyway.

Often we encounter physical parameters, whose variability is modelled by a “bell curve” … but not really: All real numbers, even unrealistic or physically impossible ones, can come out of the true \(\mathcal{N}\) distribution. This can prove most inconvenient!

  • A simple workaround, taken for several parameters in the featured model, is to follow a truncated normal distribution: samples which would fall “too far” are simply rejected and resampled. Good luck finding this in your favourite “general-purpose” model generator!



import numpy as np; import matplotlib.pyplot as plt
# set seed for reproducible figures
np.random.seed(123)#43

import scipy; from scipy import stats # for truncated normal
def truncated_normal(loc, scale, bounds, size):
    """Normal distribution truncated within bounds

    loc -- mean (“centre”) of the distribution
    scale -- standard deviation (spread or “width”) of the distribution
    bounds -- list of min and maximum
    size -- number of samples
    """
    bounds = np.array([bounds] * size)

    s = scipy.stats.truncnorm.rvs(
        (bounds[:, 0] - loc) / scale, (bounds[:, 1] - loc) / scale, loc=loc, scale=scale
    )
    return s

The point-neuron model is almost the same as <iafTauRefCell> except:

  1. it admits stimulus current as voltage;

  2. it includes a bias-current term which is different for each one neuron.

We handle the former with a <DerivedVariable> aggregating inbound V_syn, and the latter with <EdenCustomSetup>.



cell_defs = '''
<ComponentType name="VoltTauCell" extends="baseCellMembPotCap"> <!-- not really *Cap* but since baseSynapse exposes current... -->

    <Parameter name="v_reset" dimension="voltage"/>
    <Parameter name="v_thresh" dimension="voltage"/>
    <Parameter name="V_b" dimension="voltage"/>
    <Parameter name="tau_mem" dimension="time"/>
    <Parameter name="tau_refrac" dimension="time"/>

    <Dynamics>
            <StateVariable name="lastSpikeTime" dimension="time"/>
        <StateVariable name="v" dimension="voltage" exposure="v"/>
        <DerivedVariable name="V_syn" dimension="voltage" select="synapses[*]/V_syn" reduce="add" />

        <ConditionalDerivedVariable name="voltage_rate" dimension="voltage_per_time">
                    <Case condition="t .gt. lastSpikeTime + tau_refrac" value="(-v + V_syn + V_b) / tau_mem"/>
                    <Case value="0"/>
            </ConditionalDerivedVariable>

            <TimeDerivative variable="v" value="voltage_rate" />

        <OnStart> <!-- Start at resting state -->
            <StateAssignment variable="v" value="v_reset"/>
            <StateAssignment variable="lastSpikeTime" value="t - tau_refrac" />
        </OnStart>

            <OnCondition test="(v .gt. v_thresh) .and. (t .gt. lastSpikeTime + tau_refrac)">
                    <EventOut port="spike"/>
                    <StateAssignment variable="lastSpikeTime" value="t" />
                    <StateAssignment variable="v" value="v_reset" />
            </OnCondition>
    </Dynamics>
</ComponentType>

<!-- NB: V_b may be overriden with EdenCustomSetup -->
<VoltTauCell id="ExcNeuron" tau_mem="30ms" tau_refrac="3 ms" v_thresh="15 mV" v_reset="13.5 mV" V_b="0mV"  C="0.2nF" /> <!-- C is unused, see above -->
<VoltTauCell id="InhNeuron" tau_mem="30ms" tau_refrac="2 ms" v_thresh="15 mV" v_reset="13.5 mV" V_b="0mV"  C="0.2nF" />
'''

Let’s fill in the values for V_b:



rng = np.random.default_rng(seed=5)
def GetRandomIbMillivolt(rng,N,ib_mean_mvolt):
    # paper gives range of 0.05 mV but population bursts are not visible with that value
    # -> increased to 1 mV range
    return np.random.uniform(ib_mean_mvolt - 0.5, ib_mean_mvolt + 0.5, size=N)

eden_setup_lines = []
N_exc = 400; N_inh = 100

ib_vals = GetRandomIbMillivolt(rng, N=N_exc, ib_mean_mvolt=15)
# ib_vals = exc_neurons.I_b / mV # override with Brian's
eden_setup_lines += [
    f'set cell Exc all all V_b multi mV',
    'values '+' '.join([f'{v}' for v in ib_vals]) ]
ib_vals = GetRandomIbMillivolt(rng, N=N_inh, ib_mean_mvolt=15)
# ib_vals = inh_neurons.I_b / mV # override with Brian's
eden_setup_lines += [
    f'set cell Inh all all V_b multi mV',
    'values '+' '.join([f'{v}' for v in ib_vals]) ]

The synapse model is unique but not so special, other than randomised U, tau_rec, tau_facil. Until EDEN supports LEMS inheritance we’ll have to duplicate some lines:



syns_defs = '''
<ComponentType name="DynSyn" extends="baseSynapse">
    <Parameter name="A" dimension="voltage"/>
    <Parameter name="U" dimension="none"/>
    <Parameter name="tau_I" dimension="time"/>
    <Parameter name="tau_rec" dimension="time"/>
    <Constant name="AMP" dimension="current" value="1 A"/>

    <Dynamics>
        <StateVariable name="x" dimension="none"/>
        <StateVariable name="y" dimension="none"/>
        <DerivedVariable name="z" dimension="none" value="1 - x - y"/>
        <DerivedVariable name="i" exposure="i" dimension="current" value="0 * AMP" /> <!-- because baseSynapse requires it -->
        <DerivedVariable name="V_syn" exposure="V_syn" dimension="voltage" value="A*y" />

        <TimeDerivative variable="x" value="+z / tau_rec"/>
        <TimeDerivative variable="y" value="-y / tau_I"/>

        <OnStart>
            <StateAssignment variable="x" value="1"/>
            <StateAssignment variable="y" value="0"/>
        </OnStart>

        <OnEvent port="in">
            <StateAssignment variable="y" value="y + U * x"/>
            <StateAssignment variable="x" value="x - U * x"/>
        </OnEvent>

    </Dynamics>
</ComponentType>
<ComponentType name="DynSynFacil" extends="baseSynapse">
    <Parameter name="A" dimension="voltage"/>
    <Parameter name="U" dimension="none"/>
    <Parameter name="tau_I" dimension="time"/>
    <Parameter name="tau_rec" dimension="time"/>
    <Parameter name="tau_facil" dimension="time"/>
    <Constant name="AMP" dimension="current" value="1 A"/>

    <Dynamics>
        <StateVariable name="x" dimension="none"/>
        <StateVariable name="y" dimension="none"/>
        <StateVariable name="u" dimension="none"/>
        <DerivedVariable name="z" dimension="none" value="1 - x - y"/>
        <DerivedVariable name="i" exposure="i" dimension="current" value="0 * AMP" /> <!-- because baseSynapse requires it -->
        <DerivedVariable name="V_syn" exposure="V_syn" dimension="voltage" value="A*y" />

        <TimeDerivative variable="x" value="+z / tau_rec"/>
        <TimeDerivative variable="y" value="-y / tau_I"/>
        <TimeDerivative variable="u" value="-u / tau_facil"/>

        <OnStart>
            <StateAssignment variable="x" value="1"/>
            <StateAssignment variable="y" value="0"/>
            <StateAssignment variable="u" value="0"/>
        </OnStart>

        <OnEvent port="in">
            <StateAssignment variable="u" value="u + U*(1-u)"/>
            <StateAssignment variable="y" value="y + u * x"/>
            <StateAssignment variable="x" value="x - u * x"/>
        </OnEvent>
    </Dynamics>
</ComponentType>

<DynSyn      id="ToExcSyn" A="0 mV" U="0" tau_I="3 ms" tau_rec="800 ms"/>
<DynSynFacil id="ToInhSyn" A="0 mV" U="0" tau_I="3 ms" tau_rec="100 ms" tau_facil="1000 ms" />
'''

Now, let’s create all these synapses: determine pre-synaptic and post-synaptic cells, create the stochastic distributions of the various parameters, and create the set lines for <EdenCustomSetup>. (Or better, generate them while writing to the final file, to save on RAM.)



# A lot of parameters to randomise....
def GetSynsSetupDict(N_syn, tau_I, A, U, tau_rec, tau_facil=None, override=None):
    ret = {}
    ret['tau_I'] = {'vals': tau_I, 'unit': 'ms'}

    A_min = min(0.2 * A, 2 * A)
    A_max = max(0.2 * A, 2 * A)
    ret['A'] = {
        'vals': truncated_normal(
            A, 0.5 * abs(A), [A_min, A_max], size=N_syn
        ), 'unit':'mV' }
    assert not any(ret['A']['vals'] < A_min)
    assert not any(ret['A']['vals'] > A_max)
    if override: ret['A']['vals'] = override.A / mV

    U_mean, U_min, U_max = U
    ret['U'] = {'vals': truncated_normal(U_mean, 0.5 * U_mean, [U_min, U_max], size=N_syn), 'unit':''}
    assert not any(ret['U']['vals'] <= U_min)
    assert not any(ret['U']['vals'] > U_max)
    if override: ret['U']['vals'] = override.U

    tau_min = 5
    ret['tau_rec'] = {
        'vals': truncated_normal(
            tau_rec, 0.5 * tau_rec, [tau_min, np.inf], size=N_syn
        ), 'unit': 'ms' }
    assert not any(ret['tau_rec']['vals'] <= tau_min)
    if override: ret['tau_rec']['vals'] = override.tau_rec / ms

    if tau_facil:
        ret['tau_facil'] = {
            'vals': truncated_normal(
                tau_facil, 0.5 * tau_facil, [tau_min, np.inf], size=N_syn
            ), 'unit': 'ms' }
        assert not any(ret['tau_facil']['vals'] <= tau_min)
        if override: ret['tau_facil']['vals'] = override.tau_facil / ms

    return ret



# sample sparse... then randomise
def RandomUniformProjection(rng, N_from, N_to, prob, symmetric):
    M = rng.uniform(size=(N_from,N_to))
    if symmetric: M += np.eye(N_from, N_to) # exclude diagonal
    ei, ej = np.where(M < prob)
    return ei, ej

def GetSynsSetupLines(projname, property, vals, unit):
    # check if iterable https://stackoverflow.com/a/1952655
    iter8 = None
    try:
        iter8 = iter(vals)
    except TypeError: pass

    lines = [ f"set synapse {projname} all post {property} {'multi' if iter8 else vals} {unit}" ]
    if iter8: lines += [ "values "+" ".join([ str(v) for v in vals]) ]
    return lines

popsizes = {'Exc': N_exc, 'Inh': N_inh}
def get_synapses(projname, source, target, tau_I, A, U, tau_rec, tau_facil=None, override=None):
    nml_lines = [f'''\t<projection id="{projname}" presynapticPopulation="{source}" postsynapticPopulation="{target}" synapse="{'ToExcSyn' if tau_facil is None else 'ToInhSyn'}">''']
    syn_i, syn_j = RandomUniformProjection(rng, popsizes[source], popsizes[target], 0.1, (source == target))
    if override: syn_i, syn_j = override.i, override.j

    nml_lines +=  [ f'\t\t<connection id="{ii}" preCellId="{pre}" postCellId="{post}"/>' for ii, (pre, post) in enumerate(zip(syn_i, syn_j))]
    nml_lines += ['\t\t</projection>']


    setup_dict = GetSynsSetupDict(len(syn_i), tau_I, A, U, tau_rec, tau_facil, override)
    setup_lines = []
    for propname, d in setup_dict.items():
        setup_lines += GetSynsSetupLines(projname, propname, **d)
    return nml_lines, setup_lines, len(syn_i)
net_lines = []
net_ee, set_ee, n_ee = get_synapses("ee_synapses", 'Exc', 'Exc', tau_I=3, A= 1.8, U=[0.5 , 0.1  , 0.9 ], tau_rec=800                ); net_lines += net_ee; eden_setup_lines += set_ee
net_ei, set_ei, n_ei = get_synapses("ei_synapses", 'Exc', 'Inh', tau_I=3, A= 7.2, U=[0.04, 0.001, 0.07], tau_rec=100, tau_facil=1000); net_lines += net_ei; eden_setup_lines += set_ei
net_ie, set_ie, n_ie = get_synapses("ie_synapses", 'Inh', 'Exc', tau_I=3, A=-5.4, U=[0.5 , 0.1  , 0.9 ], tau_rec=800                ); net_lines += net_ie; eden_setup_lines += set_ie
net_ii, set_ii, n_ii = get_synapses("ii_synapses", 'Inh', 'Inh', tau_I=3, A=-7.2, U=[0.04, 0.001, 0.07], tau_rec=100, tau_facil=1000); net_lines += net_ii; eden_setup_lines += set_ii

Finally, let’s put it all together and run the simulation:



tabline = '\n    ' # annoyingly enough, \ may not exist at all in a f string, not even in brackets. Fixed in Python 3.12+
nml_file=f'''
<neuroml>

{cell_defs}
{syns_defs}

<!-- Here comes the network -->
<network id="Net">
    <!-- Add a population with a single cell -->
    <population id="Exc" component="ExcNeuron" size="400"/>
    <population id="Inh" component="InhNeuron" size="100"/>
    {tabline.join(net_lines)}

</network>
<Simulation id="MySim" length="5.1 s" step="100 us" target="Net">
    <EdenOutputFile id="MyOutFile" href="results.gen.txt" format="ascii_v0" sampling_interval="1 msec">
        {tabline.join([ f'<OutputColumn id="v_{i}"  quantity= "Exc[{i}]/v" output_units="mV"/>' for i in range(N_exc)])}
    </EdenOutputFile>
    <EdenOutputFile id="MyOutSyns" href="file://syns.gen.txt" format="ascii_v0" sampling_interval="1 msec">
        {tabline.join([ f'<OutputColumn id="x_{i}"  quantity= "ee_synapses[{i}]/post/x"/>' for i in range(n_ee)])}
    </EdenOutputFile>
    <EventOutputFile id="SpikesExc" fileName="spikes_exc.gen.txt" format="TIME_ID">
        {tabline.join([ f'<EventSelection id="{i}"  select="Exc[{i}]" eventPort="spike"/>' for i in range(N_exc)])}
    </EventOutputFile>
    <EventOutputFile id="SpikesInh" fileName="spikes_inh.gen.txt" format="TIME_ID">
        {tabline.join([ f'<EventSelection id="{i}"  select="Inh[{i}]" eventPort="spike"/>' for i in range(N_inh)])}
    </EventOutputFile>

    <!-- Use EdenCustomSetup ❗ -->
    <EdenCustomSetup filename="CustomSetup.txt"/>
</Simulation>
<Target component="MySim"/></neuroml>
'''
with open('Model.nml', 'wt') as f: f.write(nml_file)
setup_file = '\n'.join(eden_setup_lines)
with open('CustomSetup.txt', 'wt') as f: f.write(setup_file)

Note that we are recording the x variable of every synapse just to show their average… 400 * 400 neurons * 0.1 density * 5000 msec * 18ish ASCII bytes per number adds up to a lot! Having binary instead of ASCII output (TBD) would cut that down to 1/2 or 1/4, but complete solutions would be to either specify the average itself for recording, or to make a monstrosity with <VariableReference> powered pseudo-gap junctions into pseudo-cells for accumulation (not ideal).



import eden_simulator as eden
tra, eve = eden.runEden('Model.nml', reload_events=True, threads=1)

And sanity-check the network and see if there is synchronised firing:



v = np.array([tra[f'Exc[{i}]/v'] for i in range(N_exc)])
plt.imshow(v[:,900:], aspect='auto', interpolation='none'); plt.colorbar(label='Vm (mV)')
plt.xlabel('Time (ms)'); plt.ylabel('Neuron');
_images/paper_tsodyksuzielmarkram2000_17_0.png

And finally reproduce Figures 1a, 1c from the paper.



fig, axs = plt.subplots(2,1,sharex=True,figsize=[5,4],dpi=100) # or figsize=[11,4],dpi=200
for i in range(0,N_exc,5): # every fifth neuron for decluttering, see paper
    spikes=eve[f'Exc[{i}]']
    axs[0].plot(spikes, [i]*len(spikes), ".k", ms=1)
for i in range(0,N_inh,5):
    spikes=eve[f'Inh[{i}]']
    axs[0].plot(spikes, [i+N_exc]*len(spikes), ".k", ms=1)
axs[0].autoscale(axis='x', tight=True)
axs[0].set_ylabel("Neuron No."); axs[0].autoscale(axis='y', tight=True)
xx = np.array([tra[f'ee_synapses[{i}]/post/x'] for i in range(n_ee)])
zz = np.mean(xx,axis=0)
axs[1].plot(tra['t'],zz)
axs[1].set_ylim([0.2, 0.6]);axs[1].autoscale(axis='x', tight=True)
axs[1].set_ylabel("Rec. exc.");axs[1].set_xlabel("Time (sec)");
_images/paper_tsodyksuzielmarkram2000_19_0.png