

In [1]:

    
# import seaborn as snsf
from __future__ import print_function, division
from scipy.stats import circmean, circstd
from fakespikes.util import create_psd
from pykdf.kdf import load_kdf
from bw.util import fit_gaussian
from scipy.signal import welch
%matplotlib inline
from brian2 import *

EI populations

Increasing drive in a single population.



In [2]:

    
%run ../ie.py ie -t 3 -p 1 -q 1 --sigma .01









    



INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.15s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]






    



Starting simulation at t=0. s for a duration of 3. s
creating /var/folders/s6/8n6rytmn3xq5q0wgkxktwbs00000gn/T/scipy-type-Q4zoTj/python27_intermediate/compiler_b8fcb4d5d1597fcd5bf712a708abd317
3. s (100%) simulated in 8s



In [3]:

    
res = load_kdf('ie.hdf5')

# -
t = res['t']
dt = res['dt']
times = np.linspace(0, t, t * int(1 / float(dt)))

# -
E = res['E']
I = res['I']
lfp = res['lfp']

# -
figure(figsize=(12, 10))

subplot(311)
plot(times, lfp, label='lfp', color='grey')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

subplot(312)
plot(times, E, label='E', color='k')
plot(times, I, label='I', color='r')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

# - 
# Make and fit PSD
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=3000/2)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-3)
print(stdevs*2.355)

# -
subplot(313)
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.9)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
# plt.semilogy()
plt.tight_layout()









    



[ 2.37029708]



In [22]:

    
%run ../slidie.py slidie -t 2 --p0 0.3 --pn 1 --sigma .001

# -
res = load_kdf('slidie.hdf5')

# -
t = res['t']
dt = res['dt']
times = np.linspace(0, t, t * int(1 / float(dt)))

# -
E = res['E']
I = res['I']
lfp = res['lfp']

# -
figure(figsize=(12, 10))

subplot(311)
plot(times, lfp, label='lfp', color='grey')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
# xlim(8, 12)

subplot(312)
plot(times, E, label='E', color='k')
plot(times, I, label='I', color='r')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
# xlim(8, 12)

# - 
# Make and fit PSD
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=20000)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-3)
print(stdevs*2.355)

# -
subplot(313)
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.9)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
# plt.semilogy()
plt.tight_layout()









    



INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]






    



[ 1.79585117]

Many populations with random drives (frequencies).



In [20]:

    
%run ../mixie.py mixie -n 10 -p 2 -q 1 -s 0.5 --dt 1e-3 --seed 10 --sigma .01









    



INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.06s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]
INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]



In [21]:

    
res = load_kdf('mixie.hdf5')

# -
t = res['t']
dt = res['dt']
times = np.linspace(0, t, t * int(1 / float(dt)))

# -
E = res['E']
I = res['I']
lfp = res['lfp']

# - 
figure(figsize=(12, 10))

subplot(311)
plot(times, lfp, label='lfp', color='grey')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

subplot(312)
plot(times, E, label='E', color='k')
plot(times, I, label='I', color='r')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

# - 
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=3000)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-5)

print(stdevs*2.355)

# -
subplot(313)
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.5)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
plt.tight_layout()









    



[  1.1333496    0.94644644   0.75648667   2.09644915   0.72442748
   1.27068214   0.76593412   0.90451364   1.14702777  17.79009818
  13.22988359]

Bursting



In [19]:

    
%run ../burstie.py burstie -w 0.5 -s 1 --sigma 0.001
res = load_kdf('burstie.hdf5')

# -
t = res['t']
dt = res['dt']
times = np.linspace(0, t, t * int(1 / float(dt)))

# -
E = res['E']
I = res['I']
lfp = res['lfp']

# -
figure(figsize=(12, 10))

subplot(311)
plot(times, lfp, label='lfp', color='grey')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

subplot(312)
plot(times, E, label='E', color='k')
plot(times, I, label='I', color='r')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
xlim(0, 3)

# - 
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=3000)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-5)

print(stdevs*2.355)


# -
subplot(313)
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.5)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
plt.tight_layout()









    



INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]






    



[]

Drifting



In [17]:

    
%run ../driftie.py driftie -t 3 -d .1 --min_P 0.5 
res = load_kdf('driftie.hdf5')









    



INFO       No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.05s, trying other methods took 0.00s). [brian2.stateupdaters.base.method_choice]



In [18]:

    
# -
t = res['t']
dt = res['dt']
times = np.linspace(0, t, t * int(1 / float(dt)))

# -
E = res['E']
I = res['I']
lfp = res['lfp']

# -
figure(figsize=(12, 10))

subplot(311)
plot(times, lfp, label='lfp', color='grey')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
# xlim(0, 30)

subplot(312)
plot(times, E, label='E', color='k')
plot(times, I, label='I', color='r')
legend(loc='best')
xlabel("Time (s)")
ylabel("Rate (Hz)")
# xlim(0, 30)

# - 
# Make and fit PSD
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=1000)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-3)
# print(stdevs*2.355)


# -
subplot(313)
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.5)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
plt.tight_layout()

Coupled oscillators



In [2]:

    
%run ../kur.py kur -t 4 -n 50 -k 6 -o 30 -r 3 --sigma .1 -p 0.99 --dt 1e-3



In [3]:

    
res = load_kdf('kur.hdf5')
times = res['times']

thetas = res['thetas']
lfp = res['lfp']
waves = res['waves']

figure(figsize=(12, 3))
plot(times, thetas, color='k', alpha=0.1);



In [ ]:

    
mTheta = circmean(thetas[times > 3, :], axis=0)
r = np.ones_like(mTheta)  # Unit vectors

figure(figsize=(6, 6))
ax = plt.subplot(111, polar=True)
ax.plot(mTheta, r, '.k')
ax.plot(circmean(mTheta), 1, '.r', markersize=30, alpha=0.5)



In [ ]:

    
# Now sample the avg theta and simulate sin waves with that property, 
# each at freq range defined in the K model
# Use this to create a LFP, and PSD



In [ ]:

    
figure(figsize=(12, 3))
for n in range(res['N'])[:2]:
    wave = waves[n, :]
    if n == 0:
        plt.plot(times, wave, color='k', alpha=0.2, label='Inv. osc.')
    else:
        plt.plot(times, wave, color='k', alpha=0.2)

plt.plot(times, lfp, 'r', alpha=0.6, label='Sim. LFP')
plt.xlabel("Time (s)")
plt.ylabel("LFP (AU)")
plt.legend()
# plt.xlim(2.8, 3)
print("The oscillator frequencies are", res['omegas'])



In [ ]:

    
# - 
# Make and fit PSD
fs, psd = welch(lfp, int(1 / dt), scaling='density', nperseg=1000)

m = np.logical_and(fs > 20, fs < 40)
fs = fs[m]
psd = psd[m]
center, powers, stdevs, fit = fit_gaussian(fs, psd, 20, mph=1e-3)

print(stdevs*2.355)

# -
figure(figsize=(12, 3))
plot(fs, psd, color='k', alpha=0.9, ls='--', label='data', lw=3)
plt.plot(center, powers, 'o', alpha=0.5)
plt.plot(fs, fit, color='blue', alpha=0.3, lw=3, label='model')
plt.xlabel("Freq (Hz)")
plt.ylabel("PSD (AU)")
plt.legend()
plt.tight_layout()



In [ ]:



In [ ]: