From 8a29fe5ac7779cc08ea7af24a37fc1aa7908b390 Mon Sep 17 00:00:00 2001 From: Espen Hagen <2492641+espenhgn@users.noreply.github.com> Date: Mon, 16 Sep 2024 10:50:33 +0200 Subject: [PATCH] updated LFP plots (#229) * Accept iterable for bin sizes for pairwise spike-train correlations Fixes #216 * corrcoef plot fixes * --mpi=pmi2 not working anymore * introduce t_sim_lfp<=t_sim parameter to independently control LFP calculations * correlation analysis * Backup files before updating master * prep merge * roll back some edits * some minor plot updates * add lfp script * improved figure legibility * Update base_analysis_params.py * remove spaces * removed unused files * roll back change * add back deleted ms_figures.correlation() function * roll back edit * line change --- mesocircuit/analysis/spike_analysis.py | 4 +- mesocircuit/lfp/lfp_parameters.py | 4 +- mesocircuit/lfp/plotting.py | 39 ++++++- .../base_simulation_params.py | 4 + mesocircuit/plotting/figures.py | 25 ++--- mesocircuit/plotting/plotting.py | 6 +- mesocircuit/run/run_lfp_plotting.py | 102 ++++++++++++------ scripts/run_mesocircuit_lfps.py | 81 ++++++++++++++ 8 files changed, 207 insertions(+), 58 deletions(-) create mode 100644 scripts/run_mesocircuit_lfps.py diff --git a/mesocircuit/analysis/spike_analysis.py b/mesocircuit/analysis/spike_analysis.py index deb2c5c..426d1a7 100644 --- a/mesocircuit/analysis/spike_analysis.py +++ b/mesocircuit/analysis/spike_analysis.py @@ -951,7 +951,7 @@ def _compute_ccs_distances(X, circuit, sptrains_X, binsize_time, ccs_time_interv binsize_time Temporal resolution of sptrains_X (in ms). ccs_time_interval - Temporal bin size(s) for correlation coefficients calculation (in ms). + Temporal bin size for correlation coefficients calculation (in ms) positions_X Positions of population X. """ @@ -1018,7 +1018,7 @@ def _compute_ccs_distances(X, circuit, sptrains_X, binsize_time, ccs_time_interv extent=[circuit.net_dict['extent']] * 2, edge_wrap=True) - ccs_dic = {f'ccs_{ccs_time_interval}ms': ccs, + ccs_dic = {f'ccs_{ccs_time_interval}': ccs, 'distances_mm': distances} return ccs_dic diff --git a/mesocircuit/lfp/lfp_parameters.py b/mesocircuit/lfp/lfp_parameters.py index ee97ed8..47cd64b 100644 --- a/mesocircuit/lfp/lfp_parameters.py +++ b/mesocircuit/lfp/lfp_parameters.py @@ -212,7 +212,7 @@ def get_parameters(path_lfp_data=None, sim_dict=None, net_dict=None): savefolder=path_lfp_data, # derived parameters for CachedTopoNetwork instance network_params=dict( - simtime=sim_dict['t_presim'] + sim_dict['t_sim'], + simtime=sim_dict['t_presim'] + sim_dict['t_sim_lfp'], dt=2**-2, # sim_dict['sim_resolution'], # run at 4kHz; we're downsampling # outputs anyway. @@ -288,7 +288,7 @@ def get_parameters(path_lfp_data=None, sim_dict=None, net_dict=None): 'lambda_f': 100, 'dt': 2**-2, # sim_dict['sim_resolution'], 'tstart': 0, - 'tstop': sim_dict['t_presim'] + sim_dict['t_sim'], + 'tstop': sim_dict['t_presim'] + sim_dict['t_sim_lfp'], 'verbose': False, } )) diff --git a/mesocircuit/lfp/plotting.py b/mesocircuit/lfp/plotting.py index 14af5d7..7cbd948 100644 --- a/mesocircuit/lfp/plotting.py +++ b/mesocircuit/lfp/plotting.py @@ -462,7 +462,8 @@ def layout_illustration( def plot_single_channel_lfp_data( ax, PS, fname, title='LFP', ylabel=r'$\Phi$ (mV)', T=[ 500, 550], CONTACTPOS=( - (-200, 200), (200, -200)), subtract_mean=True): + (-200, 200), (200, -200)), subtract_mean=True, + report_corrcoefs=False): ''' Arguments --------- @@ -476,6 +477,8 @@ def plot_single_channel_lfp_data( x and y coordinate of electrode contact points in PS.electrodeParams subtract_mean: bool if True, subtract mean value trace + report_corrcoefs: bool + if True, print pairwise correlations in axes title ''' with h5py.File(fname, 'r') as f: srate = f['srate'][()] @@ -491,6 +494,20 @@ def plot_single_channel_lfp_data( ax.plot(tvec, data[tinds] - data[tinds].mean(), f'C{i}') else: ax.plot(tvec, data[tinds], f'C{i}') + + if report_corrcoefs: + with h5py.File(fname, 'r') as f: + tind = int(T[0] * srate / 1000) + CONTACTS = [] + for i, CPOS in enumerate(CONTACTPOS): + CONTACT = (PS.electrodeParams['x'] == CPOS[0] + ) & (PS.electrodeParams['y'] == CPOS[1]) + CONTACTS += np.where(CONTACT) + CONTACTS = np.array(CONTACTS).ravel() + data = f['data'][()][CONTACTS, tind:] + cc = np.corrcoef(data) + print(f'{title} correlations', CONTACTS + 1, CONTACTPOS, cc) + #else: ax.set_title(title) ax.set_ylabel(ylabel) remove_axis_junk(ax) @@ -501,7 +518,8 @@ def plot_single_channel_csd_data( title='CSD', ylabel=r'CSD ($\frac{\mathrm{nA}}{\mathrm{µm}^3}$)', T=[500, 550], - CONTACTPOS=((-200, 200), (200, -200))): + CONTACTPOS=((-200, 200), (200, -200)), + report_corrcoefs=False): # CSD bin edges X, Y, Z = np.meshgrid(PS.CSDParams['x'], PS.CSDParams['y'], @@ -523,6 +541,16 @@ def plot_single_channel_csd_data( CONTACT = (Xmid == CPOS[0]) & (Ymid == CPOS[1]) data = f['data'][()][CONTACT, ].flatten() ax.plot(tvec, data[tinds] - data[tinds].mean(), f'C{i}') + if report_corrcoefs: + datas = [] + with h5py.File(fname, 'r') as f: + tind = int(T[0] * srate / 1000) + for i, CPOS in enumerate(CONTACTPOS): + CONTACT = (Xmid == CPOS[0]) & (Ymid == CPOS[1]) + datas += [f['data'][()][CONTACT, ].flatten()[tind:]] + cc = np.corrcoef(np.array(datas)) + print(f'{title} correlations', cc) + ax.set_title(title) ax.set_ylabel(ylabel) remove_axis_junk(ax) @@ -1244,7 +1272,8 @@ def plot_coherence_vs_distance_vs_frequency( method='mlab', phase_coherence=False, title='LFP', - show_cbar=True + show_cbar=True, + clim=[0, 0.25] ): with h5py.File(fname, 'r') as f: @@ -1279,8 +1308,8 @@ def plot_coherence_vs_distance_vs_frequency( unique * 1E-3, # [µm] -> [mm] unit conversion chfreqs, means, - vmin=0, - vmax=0.25, + vmin=clim[0], + vmax=clim[1], cmap='viridis', shading='auto') diff --git a/mesocircuit/parameterization/base_simulation_params.py b/mesocircuit/parameterization/base_simulation_params.py index 7cbe8c4..eea6937 100644 --- a/mesocircuit/parameterization/base_simulation_params.py +++ b/mesocircuit/parameterization/base_simulation_params.py @@ -30,6 +30,10 @@ # resolution of the simulation (in ms) 'sim_resolution': 0.1, + # simulation time for LFP predictions (in ms), + # should be equal to or shorter than 't_sim' above. + 't_sim_lfp': 1000.0, + # list of recording devices, default is 'spike_recorder'. A 'voltmeter' can # be added to record membrane voltages of the neurons. Nothing will be # recorded if an empty list is given. diff --git a/mesocircuit/plotting/figures.py b/mesocircuit/plotting/figures.py index 1053a44..1d431da 100644 --- a/mesocircuit/plotting/figures.py +++ b/mesocircuit/plotting/figures.py @@ -4,14 +4,12 @@ Definition of default figures plotted with plotting.py. """ -import mesocircuit.plotting.plotting as plot +from mesocircuit.plotting import plotting as plot from mesocircuit.parameterization import helpers_analysis as helpana import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt import numpy as np import h5py -import matplotlib -matplotlib.use('Agg') def parameters(circuit): @@ -313,14 +311,11 @@ def statistics_overview(circuit, all_FRs, all_LVs, all_CCs_distances, all_PSDs): for X in all_CCs_distances: if isinstance(all_CCs_distances[X], h5py._hl.group.Group): try: - iter(circuit.ana_dict['ccs_time_interval']) - ccs_time_interval = circuit.ana_dict['ccs_time_interval'][0] - all_CCs[X] = all_CCs_distances[X][f'ccs_{ccs_time_interval}ms'] - - print('CCs in statistics_overview use only first value in list of time intervals: ' + - f'{ccs_time_interval}') - except TypeError: - all_CCs[X] = all_CCs_distances[X][f'ccs_{circuit.ana_dict["ccs_time_interval"]}ms'] + for i, ccs_time_interval in enumerate(iter(circuit.ana_dict['ccs_time_interval'])): + all_CCs[X] = all_CCs_distances[X][f'ccs_{ccs_time_interval}'] + break # plot statistics overview using only first value in list for now + except TypeError as _: + all_CCs[X] = all_CCs_distances[X][f'ccs_{circuit.ana_dict["ccs_time_interval"]}'] else: all_CCs[X] = np.array([]) @@ -365,20 +360,20 @@ def _corrcoef_distance(key_ccs): xlabel='distance (mm)', ylabel=r'$CC$') return fig - + try: - # one figure per time interval if a list is provided + # stats plot for different time bins for i, ccs_time_interval in enumerate(iter(circuit.ana_dict['ccs_time_interval'])): key_ccs = f'ccs_{ccs_time_interval}' fig = _corrcoef_distance(key_ccs) plot.savefig(circuit.data_dir_circuit, circuit.plot_dict['extension'], - f'corrcoef_distance_{key_ccs}ms') + f'corrcoef_distance_{key_ccs}ms') except TypeError as _: ccs_time_interval = circuit.ana_dict['ccs_time_interval'] key_ccs = f'ccs_{ccs_time_interval}' fig = _corrcoef_distance(key_ccs) plot.savefig(circuit.data_dir_circuit, circuit.plot_dict['extension'], - f'corrcoef_distance_{key_ccs}ms') + f'corrcoef_distance_{key_ccs}ms') return diff --git a/mesocircuit/plotting/plotting.py b/mesocircuit/plotting/plotting.py index d19c627..e21697c 100644 --- a/mesocircuit/plotting/plotting.py +++ b/mesocircuit/plotting/plotting.py @@ -23,7 +23,7 @@ import matplotlib from scipy.optimize import curve_fit from mesocircuit.helpers.io import load_h5_to_sparse_X -matplotlib.use('Agg') + matplotlib.rcParams.update(rcParams) # initialize MPI @@ -1368,7 +1368,7 @@ def plotfunc_CCs_distance( return distances = data[X]['distances_mm'][:max_num_pairs] - ccs = data[X][key_ccs + 'ms'][:max_num_pairs] + ccs = data[X][key_ccs][:max_num_pairs] # loop for reducing zorder-bias blocksize = int(len(distances) / nblocks) @@ -1517,7 +1517,7 @@ def savefig( data_dir_circuit Path to data directory extension - File extension. + file extension filename Name of the file. eps_conv diff --git a/mesocircuit/run/run_lfp_plotting.py b/mesocircuit/run/run_lfp_plotting.py index c90f397..94dbb0a 100644 --- a/mesocircuit/run/run_lfp_plotting.py +++ b/mesocircuit/run/run_lfp_plotting.py @@ -76,14 +76,12 @@ 'Population', 'run', 'collect']) ], ignore_index=True) -df.to_csv('simstats.csv', index=False) +df.to_csv(os.path.join(path_lfp_data, 'simstats.csv'), index=False) - -df = pd.read_csv('simstats.csv') # sum df['sum'] = np.round( df[['CachedNetwork', 'Population', 'run', 'collect'] - ].sum(axis=1)).astype(int) + ].astype(float).sum(axis=1)).astype(int) # per second of total simulation duration df['per_s'] = np.round(df['sum'] / (sim_dict['t_presim'] + @@ -101,9 +99,6 @@ ax.set_ylabel('time (s)') ax.set_title('simulation time') -print(df[['y', 'per_s']]) -print(df[['per_s']].to_numpy().flatten().tolist()) - ######################################################################### # Create an object representation containing the spiking activity of the @@ -156,20 +151,23 @@ T = [sim_dict['t_presim'], sim_dict['t_presim'] + 500] fname = os.path.join(path_lfp_data, PS.electrodeFile) lfpplt.plot_single_channel_lfp_data(axes[0], PS, fname, - T=T, CONTACTPOS=CONTACTPOS) + T=T, CONTACTPOS=CONTACTPOS, + report_corrcoefs=True) plt.setp(axes[0].get_xticklabels(), visible=False) # Figure 7C: plot CSD in same channel fname = os.path.join(path_lfp_data, PS.CSDFile) lfpplt.plot_single_channel_csd_data(axes[1], PS, fname, - T=T, CONTACTPOS=CONTACTPOS) + T=T, CONTACTPOS=CONTACTPOS, + report_corrcoefs=True) plt.setp(axes[1].get_xticklabels(), visible=False) # Figure 7D: plot MUA in same channel fname = os.path.join(path_lfp_data, PS.MUAFile) lfpplt.plot_single_channel_lfp_data(axes[2], PS, fname, T=T, CONTACTPOS=CONTACTPOS, - title='MUA', ylabel=r'$s^{-1}$') + title='MUA', ylabel=r'$s^{-1}$', + report_corrcoefs=True) axes[2].set_xlabel('time (ms)') fig.savefig(os.path.join(path_fig_files, 'signal_timeseries_I.pdf')) @@ -181,7 +179,9 @@ plot_dict['fig_width_2col'], plot_dict['fig_width_1col'])) axes = [] -gs = GridSpec(3, 4, wspace=0.7, hspace=0.3) +gs = GridSpec(3, 4, + left=0.04, right=0.975, bottom=0.08, top=0.925, + wspace=0.6, hspace=0.25) # Figure 7A CONTACTPOS = ((600, 600), (600, 1000), (-1400, -1400)) @@ -341,7 +341,8 @@ figsize=(plot_dict['fig_width_2col'], plot_dict['fig_width_1col'] * 1.5), sharex=True) -fig.subplots_adjust(wspace=0.15, hspace=0.15) +fig.subplots_adjust(left=0.08, right=0.975, bottom=0.05, top=0.95, + wspace=0.15, hspace=0.15) axes = axes.flatten() # Fig 8 A spike train correlations @@ -364,7 +365,7 @@ fname = os.path.join( circuit.data_dir, 'upscaled_CCs_only', - '3be54d189b4f3a23c4859e0aea6b5fa8', + 'b09dcfc3953647c1a29c007fe8dde5fa', 'processed_data', 'all_CCs_distances.h5') assert os.path.isfile(fname) @@ -377,9 +378,10 @@ f'for figure_08.pdf in file {fname}') with h5py.File(fname, 'r') as f: - for i, (X, n_pairs) in enumerate(zip(['L23E', 'L23I'], [40, 10])): + for i, (X, n_pairs) in enumerate(zip(['L23E', 'L23I'], [512, 512])): plot.plotfunc_CCs_distance( ax=ax, X=X, i=i, data=f, + key_ccs='ccs_2.0ms', pop_colors=plot_dict['pop_colors'], max_num_pairs=n_pairs, markersize_scale=0.4, @@ -388,7 +390,7 @@ # add entries with NaNs to mask for i, X in enumerate(['L23E', 'L23I']): - c = f[X]['ccs'][()] + c = f[X]['ccs_2.0ms'][()] mask = c != np.nan bins = np.linspace(c[mask].mean() - c[mask].std() * 1.5, @@ -448,10 +450,10 @@ lfpplt.plot_coherence_vs_frequency( ax, PS, fname, title, colors=plt.get_cmap('viridis', 5), - NFFT=256, - noverlap=196, + NFFT=ana_dict['psd_NFFT'], + noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), method='mlab', - tbin=0.5, + tbin=ana_dict['binsize_time'], TRANSIENT=sim_dict['t_presim']) fig.savefig(os.path.join(path_fig_files, 'signal_coherence.pdf')) @@ -465,8 +467,8 @@ lfpplt.plot_coherence_vs_distance( ax, PS, fname, max_inds=np.array([2, 6, 16, 26, 38]), - NFFT=256, noverlap=196, - method='mlab', tbin=0.5, + NFFT=ana_dict['psd_NFFT'], noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), + method='mlab', tbin=ana_dict['binsize_time'], fit_exp=fit_exp, phase_coherence=False) if fit_exp: @@ -488,10 +490,10 @@ for ax, fname, title in zip(axes, fnames, titles): lfpplt.plot_coherence_vs_distance_vs_frequency( fig, ax, PS, fname, - NFFT=256, - noverlap=196, + NFFT=ana_dict['psd_NFFT'], + noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), method='mlab', - tbin=0.5, + tbin=ana_dict['binsize_time'], TRANSIENT=sim_dict['t_presim'], title=title ) @@ -505,40 +507,78 @@ fig, axes = plt.subplots(2, 3, figsize=(plot_dict['fig_width_2col'], plot_dict['fig_width_1col'] * 1.5), sharex='row', sharey='row') -fig.subplots_adjust(hspace=0.4, wspace=0.4) +fig.subplots_adjust(left=0.065, right=0.925, bottom=0.05, top=0.95, + hspace=0.3, wspace=0.3) fnames = [os.path.join(path_lfp_data, PS.electrodeFile), os.path.join(path_lfp_data, PS.CSDFile), os.path.join(path_lfp_data, PS.MUAFile)] +def _get_clims(): + # assess the + with h5py.File(fnames[0], 'r') as f: + Fs = f['srate'][()] + T0 = int(Fs * sim_dict['t_presim'] / 1000) # t < T0 transient + shape = f['data'].shape + if len(shape) > 2: + # flatten all but last axis + data = f['data'][()].reshape((-1, shape[-1]))[:, T0:] + else: + data = f['data'][()][:, T0:] + + data_x = data + data_y = data + + r, c, chfreqs, mask = lfpplt.get_data_coherence( + data_x=data_x, data_y=data_y, srate=Fs, + positions_x=PS.electrodeParams['x'], + positions_y=PS.electrodeParams['y'], + tbin=ana_dict['binsize_time'], NFFT=ana_dict['psd_NFFT'], noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), + method='mlab', + phase_coherence=False) + + unique = np.unique(r[mask]) + means = [] + for d in unique: + if d == 0: + continue + means += [c[mask][r[mask] == d].mean(axis=0)] + means = np.array(means).T + + return [0, means[chfreqs >= 10, :].max()] + +clim = _get_clims() + # panels A-C titles = ['LFP', 'CSD', 'MUA'] for i, (ax, fname, title) in enumerate(zip(axes[0, ], fnames, titles)): lfpplt.plot_coherence_vs_frequency( ax, PS, fname, title, colors=plt.get_cmap('viridis', 5), - NFFT=256, - noverlap=196, + NFFT=ana_dict['psd_NFFT'], + noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), method='mlab', - tbin=0.5, + tbin=ana_dict['binsize_time'], TRANSIENT=sim_dict['t_presim'], show_legend=True if i == 0 else False) if i > 0: plt.setp(ax.get_yticklabels(), visible=False) ax.set_ylabel('') ax.axis('tight') +ax.set_ylim(clim) # panels D-F for i, (ax, fname, title) in enumerate(zip(axes[1, :], fnames, titles)): lfpplt.plot_coherence_vs_distance_vs_frequency( fig, ax, PS, fname, - NFFT=256, - noverlap=196, + NFFT=ana_dict['psd_NFFT'], + noverlap=int(ana_dict['psd_NFFT'] * 3 // 4), method='mlab', - tbin=0.5, + tbin=ana_dict['binsize_time'], TRANSIENT=sim_dict['t_presim'], title=title, - show_cbar=(i == 2) + show_cbar=(i == 2), + clim=[clim[0], clim[1] * 0.5] ) if i > 0: plt.setp(ax.get_yticklabels(), visible=False) diff --git a/scripts/run_mesocircuit_lfps.py b/scripts/run_mesocircuit_lfps.py new file mode 100644 index 0000000..fbc935e --- /dev/null +++ b/scripts/run_mesocircuit_lfps.py @@ -0,0 +1,81 @@ +"""PyNEST Mesocircuit: Run Mesocircuit +-------------------------------------- + +Main example script to run a simulation of the mesocircuit with NEST and +subsequently analyze and plot the results. +The simulation of the spiking neural network can be followed by an LFP +simulation and the corresponding postprocessing and plotting. +""" + +############################################################################### +# The script `mesocircuit_framework` contains the main functionality for +# parameter evalutation and job execution. +# The mesocircuit framework allows for the evaluation of parameter spaces, but +# here we only simulate one individual parameter combination. +# Several interesting parameter combinations (overwriting default values) are +# collected in a dictionary in the script `parametersets` for convenience. + +from mesocircuit import mesocircuit_framework as mesoframe +import parametersets + +################################################################################ +# Here, we choose the parameter set `mesocircuit_MAMV1` which is the default +# model. It is an upscaled version of the microcircuit representing area V1 of +# the Multi-Area Model (Schmidt and van Albada, 2018). +# 'mesocircuit_MAMV1_evoked` applies a thalamic stimulus to the center of the +# same model. +# For local testing, `local_microcircuit_PD` and `local_mesocircuit_PD` are good +# choices. These models base on the original microcircuit +# (Potjans and Diesmann, 2014) and are downscaled for execution on a laptop. + +name = 'mesocircuit_MAMV1' +# name = 'local_mesocircuit_PD' +custom_params = parametersets.ps_dicts[name] + +# set long dur +custom_params['sim_dict'] = dict(t_sim=10000, t_sim_lfp=5000) +custom_params['sys_dict']['hpc'].update( + {'analysis_and_plotting': {'wall_clock_time': '02:00:00'}}) + +################################################################################ +# Next, we instantiate a `MesocircuitExperiment` with the custom parameters. +# The argument `name` can be chosen freely; here we just use the name of the +# parameter set. +# Upon instantiation, data directories are created, derived parameters +# calculated, and job scripts written. +# If an already existing experiment should be loaded, a class can be +# instantiated with the arguments of the existing `name` and `load=True`. + +meso_exp = mesoframe.MesocircuitExperiment(name, custom_params) + + +################################################################################ +# A `MesocircuitExperiment` provides an overview over all the parameter +# combinations it is holding (`parameterview`) and a list of all the individual +# model instances of class `Mesocircuit`` (`circuits`). + +print(meso_exp.parameterview) +print(meso_exp.circuits) + +################################################################################ +# For each Mesocircuit jobs can finally be launched. +# All provided jobs are run one after the other. `analysis_and_plotting` is for +# convenience combined into one job, but `analysis` and `plotting` can also be +# handled as individual jobs. +# For running the full model on an HPC cluster, `machine='hpc'` is required for +# submitting batch scripts via slurm; for a local test run `machine='local'` +# should be selected. + +for circuit in meso_exp.circuits: + circuit.run_jobs( + jobs=[ + # 'network', + # 'analysis', + # 'plotting', + # 'lfp_simulation', + # 'lfp_postprocess', + 'lfp_plotting', + ], + machine='hpc', + # machine='local' + )