Combining alpha/beta with an ERP

[tbardouille]

Hi there,

I'm interested in looking at the effect of an ERP on the alpha/beta complexes in HNN simulations. I've started by simulating a dipole for a network that combines (1) ongoing alpha/beta (as per the Example 2) for 1000 ms, with (2) an ERP (as per Example 1) beginning at 400 ms.

What I've found is that the ERP signal is much larger (about 1 order of magnitude) than the alpha/beta signals. This makes is difficult for me to further investigate here.

I'd like to adjust the network so that the two sets of signals (ERP and alpha/beta) are more closely matched in magnitude (as per experimental findings), while maintaining the same drive characteristics in the model. I'm not sure how to go about doing this. It seems like I need a smaller ERP network inside of a larger alpha/beta network, but I don't know how to start on that.

I've attached my development code (additional post) and an image of the output (attached). Do you have any advice for me?

Thanks,
Tim.

[tbardouille]

import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, law_2021_model
from hnn_core.viz import plot_dipole
import matplotlib.pyplot as plt

# Generate alpha and beta bursts with an ERP co-occurring based on examples 1 & 2 at hnn_core

# Authors: Mainak Jas <mjas@mgh.harvard.edu>
#          Sam Neymotin <samnemo@gmail.com>
#          Nick Tolley <nicholas_tolley@brown.edu>
#          Christopher Bailey <bailey.cj@gmail.com>
#          Blake Caldwell <blake_caldwell@brown.edu>
#          Tim Bardouille <tim.bardouille@dal.ca>

trialDuration = 1000.
erpDelay = 400.

# Instantiate network
net = jones_2009_model()
#net = law_2021_model()

## Add 10 Hz Bursty Proximal Drive
location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 5.4e-5, 'L5_pyramidal': 5.4e-5}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net.add_bursty_drive(
    'alpha_prox', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,
    spike_isi=10, n_drive_cells=10, location=location,
    weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=284)

## 10 Hz Bursty Distal Drive
location = 'distal'
burst_std = 15
weights_ampa_d = {'L2_pyramidal': 5.4e-5, 'L5_pyramidal': 5.4e-5}
syn_delays_d = {'L2_pyramidal': 5., 'L5_pyramidal': 5.}
net.add_bursty_drive(
    'alpha_dist', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,
    spike_isi=10, n_drive_cells=10, location=location,
    weights_ampa=weights_ampa_d, synaptic_delays=syn_delays_d, event_seed=296)

# Distal evoked drive
weights_ampa_d1 = {'L2_basket': 0.006562, 'L2_pyramidal': .000007,
                   'L5_pyramidal': 0.142300}
weights_nmda_d1 = {'L2_basket': 0.019482, 'L2_pyramidal': 0.004317,
                   'L5_pyramidal': 0.080074}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                      'L5_pyramidal': 0.1}
net.add_evoked_drive(
    'evdist1', mu=erpDelay+63.53, sigma=3.85, numspikes=1, weights_ampa=weights_ampa_d1,
    weights_nmda=weights_nmda_d1, location='distal',
    synaptic_delays=synaptic_delays_d1, event_seed=274)

# Proximal drive 1
weights_ampa_p1 = {'L2_basket': 0.08831, 'L2_pyramidal': 0.01525,
                   'L5_basket': 0.19934, 'L5_pyramidal': 0.00865}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                        'L5_basket': 1., 'L5_pyramidal': 1.}
# all NMDA weights are zero; pass None explicitly
net.add_evoked_drive(
    'evprox1', mu=erpDelay+26.61, sigma=2.47, numspikes=1, weights_ampa=weights_ampa_p1,
    weights_nmda=None, location='proximal',
    synaptic_delays=synaptic_delays_prox, event_seed=544)

# Proximal drive 2. NB: only AMPA weights differ from first
weights_ampa_p2 = {'L2_basket': 0.000003, 'L2_pyramidal': 1.438840,
                   'L5_basket': 0.008958, 'L5_pyramidal': 0.684013}
# all NMDA weights are zero; omit weights_nmda (defaults to None)
net.add_evoked_drive(
    'evprox2', mu=erpDelay+137.12, sigma=8.33, numspikes=1,
    weights_ampa=weights_ampa_p2, location='proximal',
    synaptic_delays=synaptic_delays_prox, event_seed=814)

# Simulate dipole timecourse
dpl = simulate_dipole(net, tstop=trialDuration, n_trials=1)

fig, ax = plt.subplots(1,1,figsize=(8,6))
hnn_core.viz.plot_dipole(dpl, ax=ax, decim=40)
fig.savefig('alphaBeta_and_ERP_fig1.png')
plt.close()

dplLP = dpl[0].copy()
dplLP.savgol_filter(30)

fig, ax = plt.subplots(1,1,figsize=(8,6))
hnn_core.viz.plot_dipole(dplLP, ax=ax, decim=40)
fig.savefig('alphaBeta_and_ERP_fig2.png')
plt.close()

[jasmainak]

@tbardouille perhaps this is a useful starting point:

https://jonescompneurolab.github.io/hnn-core/stable/auto_examples/workflows/plot_simulate_beta.html#sphx-glr-auto-examples-workflows-plot-simulate-beta-py

?

[tbardouille]

Hi Jas,

Thanks for replying. Using the example 5 code (Law), the ERP signal (row 2 below) is still about 1 order of magnitude than the beta signal (see row 1 below). This doesn't look like a way to address what I'm trying to do.

Is there a way to make the drives that generate the ERP only hit a subset of the network, thus creating a smaller magnitude response? Or maybe you have some other thoughts on this?

Best,
Tim.

[tbardouille]

My code for the above:

# Authors: 
#   Nick Tolley <nicholas_tolley@brown.edu>
#   Tim Bardouille <tim.bardouille@dal.ca>

from hnn_core import simulate_dipole, law_2021_model, jones_2009_model
from hnn_core.viz import plot_dipole
import hnn_core
import matplotlib.pyplot as plt
import numpy as np

net = law_2021_model()

def add_erp_drives(net, stimulus_start):
    # Distal evoked drive
    weights_ampa_d1 = {'L2_basket': 0.0005, 'L2_pyramidal': 0.004,
                       'L5_pyramidal': 0.0005}
    weights_nmda_d1 = {'L2_basket': 0.0005, 'L2_pyramidal': 0.004,
                       'L5_pyramidal': 0.0005}
    syn_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                     'L5_pyramidal': 0.1}
    net.add_evoked_drive(
        'evdist1', mu=70.0 + stimulus_start, sigma=0.0, numspikes=1,
        weights_ampa=weights_ampa_d1, weights_nmda=weights_nmda_d1,
        location='distal', synaptic_delays=syn_delays_d1, event_seed=274)

    # Two proximal drives
    weights_ampa_p1 = {'L2_basket': 0.002, 'L2_pyramidal': 0.0011,
                       'L5_basket': 0.001, 'L5_pyramidal': 0.001}
    syn_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                       'L5_basket': 1., 'L5_pyramidal': 1.}

    # all NMDA weights are zero; pass None explicitly
    net.add_evoked_drive(
        'evprox1', mu=25.0 + stimulus_start, sigma=0.0, numspikes=1,
        weights_ampa=weights_ampa_p1, weights_nmda=None,
        location='proximal', synaptic_delays=syn_delays_prox, event_seed=544)

    # Second proximal evoked drive. NB: only AMPA weights differ from first
    weights_ampa_p2 = {'L2_basket': 0.005, 'L2_pyramidal': 0.005,
                       'L5_basket': 0.01, 'L5_pyramidal': 0.01}
    # all NMDA weights are zero; omit weights_nmda (defaults to None)
    net.add_evoked_drive(
        'evprox2', mu=135.0 + stimulus_start, sigma=0.0, numspikes=1,
        weights_ampa=weights_ampa_p2, location='proximal',
        synaptic_delays=syn_delays_prox, event_seed=814)

    return net

def add_beta_drives(net, beta_start):
    # Distal Drive
    weights_ampa_d1 = {'L2_basket': 0.00032, 'L2_pyramidal': 0.00008,
                       'L5_pyramidal': 0.00004}
    syn_delays_d1 = {'L2_basket': 0.5, 'L2_pyramidal': 0.5,
                     'L5_pyramidal': 0.5}
    net.add_bursty_drive(
        'beta_dist', tstart=beta_start, tstart_std=0., tstop=beta_start + 50.,
        burst_rate=1., burst_std=10., numspikes=2, spike_isi=10,
        n_drive_cells=10, location='distal', weights_ampa=weights_ampa_d1,
        synaptic_delays=syn_delays_d1, event_seed=290)

    # Proximal Drive
    weights_ampa_p1 = {'L2_basket': 0.00004, 'L2_pyramidal': 0.00002,
                       'L5_basket': 0.00002, 'L5_pyramidal': 0.00002}
    syn_delays_p1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,
                     'L5_basket': 1.0, 'L5_pyramidal': 1.0}

    net.add_bursty_drive(
        'beta_prox', tstart=beta_start, tstart_std=0., tstop=beta_start + 50.,
        burst_rate=1., burst_std=20., numspikes=2, spike_isi=10,
        n_drive_cells=10, location='proximal', weights_ampa=weights_ampa_p1,
        synaptic_delays=syn_delays_p1, event_seed=300)
    return net

beta_start, stimulus_start = 50.0, 125.0
net_beta = net.copy()
net_beta = add_beta_drives(net_beta, beta_start)

net_beta_erp = net_beta.copy()
net_beta_erp = add_erp_drives(net_beta_erp, stimulus_start)

net_erp = net.copy()
net_erp = add_erp_drives(net_erp, stimulus_start)

dpls_beta = simulate_dipole(net_beta, tstop=400)
dpls_erp = simulate_dipole(net_erp, tstop=400)
dpls_beta_erp = simulate_dipole(net_beta_erp, tstop=400)

fig, axes = plt.subplots(4, 1, sharex=True, figsize=(7, 7),
                         constrained_layout=True)
net_beta.cell_response.plot_spikes_hist(ax=axes[0], show=False)
axes[0].set_title('Beta Event Generation')
plot_dipole(dpls_beta, ax=axes[1], layer='agg', tmin=1.0, color='b', show=False)
net_beta.cell_response.plot_spikes_raster(ax=axes[2], show=False)
axes[2].set_title('Spike Raster')

# Create a fixed-step tiling of frequencies from 1 to 40 Hz in steps of 1 Hz
freqs = np.arange(10., 60., 1.)
dpls_beta[0].plot_tfr_morlet(freqs, n_cycles=7, ax=axes[3])
fig.savefig('lawExample_fig0.png')
plt.close()


fig, ax = plt.subplots(3,1,figsize=(8,6), sharex=True)
hnn_core.viz.plot_dipole(dpls_beta, ax=ax[0], decim=40)
ax[0].legend(['Beta Only'])
hnn_core.viz.plot_dipole(dpls_erp, ax=ax[1], decim=40)
ax[1].legend(['ERP Only'])
hnn_core.viz.plot_dipole(dpls_beta_erp, ax=ax[2], decim=40)
ax[2].legend(['Beta and ERP'])
fig.savefig('lawExample_fig1.png')
plt.close()

[jasmainak]

Is there a way to make the drives that generate the ERP only hit a subset of the network, thus creating a smaller magnitude response? Or maybe you have some other thoughts on this?

Have you tried setting the probability parameter in add_bursty_drive? Another possibility is to delete the connection from net.connectivity corresponding to the drive using pick_connectivity and adding back the desired connectivity pattern using add_connection.

I guess you could also play around with the synaptic weights but I'm not sure if it's the same thing ... @rythorpe may have some further thoughts, he has better intuitions on how the network behaves than me :)

[rythorpe]

Hey @tbardouille, this is a really interesting problem - one that we haven't quite resolved in our own heads yet because the Sherman et al model for beta events assumes subthreshold postsynaptic primary currents for almost the entire network. As I think you're hinting at, it's unknown whether such beta events emerge from a subpopulation within a neocortical area that specifically doesn't respond much to evoked activation (thus avoiding the confound of pronounced network spiking), or if the coincident prox+dist drives we hypothesize produce subthreshold beta events are only part of the story.

All that said, I like your idea about targeting a subset of the network with evoked drives in order to reduce the relative balance of evoked response to beta event. The solution I'd recommend might open a can of worms, but it involves resetting the space_constant parameter for both the local network connectivity and the evoked drive inputs. The goal would be to recreate the sequence of spiking that propagates throughout the network in response to the evoked drives, but just in a smaller subset of cells. The complicated part is that you might need to adjust the space_constant (sometimes referred to as lamtha in the code, e.g., connection['nc_dict']['lamtha']) and synaptic weight for each type of local network and evoked drive connection.

One other issue you might run into if you try implementing my suggestion above is that currently, both the bursty and evoked drives target the center of the network, thus exciting the same subpopulation. This might limit your ability to observe a robust beta event coincident with a robust evoked response. There might be a way to hack this, but we should probably meet offline to discuss. LMK if you'd like to chat more - I'm happy to setup a zoom call or something.

[tbardouille]

Hey @jasmainak and @rythorpe,

Yeah, it would be helpful to chat more about this. I'm just wrapping up teaching duties, and then will try some of what you've suggested with Lindsey's help. So, you can expect to hear from us in early May with our findings. It sounds like the solution is not obvious/simple. In a way, that is good for us as there is a real chance for collaboration and useful output for the research community! :)

Can you please add @lindseypower to this discussion?

Thanks,
Tim.

[rythorpe]

Sounds good @tbardouille and @lindseypower - let's chat more about this soon!