Abstract¶

print('-------> Tested under the following settings:')
import scipy
print('SciPy version :', scipy.__version__)
import numpy
print('NumPy version :', numpy.__version__)
import matplotlib
print('Matplotlib version :', matplotlib.__version__)
import brian2
print('Brian2 version :', brian2.__version__)

# displaying figures inline:
%matplotlib inline 

# the parameters are written within a library of neuron parameters
from modeling_mesoscopic_dynamics.single_cell_models.cell_library import get_neuron_params
RS_params = get_neuron_params('RS-cell', name='RS-cell', SI_units=False)
print(RS_params)

# plotting the responses to current steps
from modeling_mesoscopic_dynamics.single_cell_models.step_response import make_model_figure
make_model_figure('RS-cell', I0=200e-12);

# Now making the numerical simulations at various levels of both excitatory and inhibitory inputs
# note the rescaling of the excitatory input to insure that the response remains in the <30 range: it focuses on low range when inhibition is low and on a larger range when inhibition is larger
%run  modeling_mesoscopic_dynamics/transfer_functions/tf_simulation.py RS-cell CONFIG1 -s --SEED 4 --tstop 10 --discret_Fe 10  --discret_Fi 10 --max_Fe 30
# paper's value -> 5-10 minutes sim

from modeling_mesoscopic_dynamics.transfer_functions.theoretical_tools import make_fit_from_data
P_RS = make_fit_from_data('data/RS-cell_CONFIG1.npy', with_square_terms=True, verbose=False)

# now displaying the numerical values and the analytical function
from modeling_mesoscopic_dynamics.transfer_functions.plots import make_exc_inh_fig
make_exc_inh_fig('data/RS-cell_CONFIG1.npy', P=P_RS)

# the parameters are written within a library of neuron parameters
from modeling_mesoscopic_dynamics.single_cell_models.cell_library import get_neuron_params
FS_params = get_neuron_params('FS-cell', name='FS-cell', SI_units=False)
print(FS_params)

# plotting the responses to current steps
from modeling_mesoscopic_dynamics.single_cell_models.step_response import make_model_figure
make_model_figure('FS-cell', I0=200e-12);

# Now making the numerical simulations at various levels of both excitatory and inhibitory inputs
%run modeling_mesoscopic_dynamics/transfer_functions/tf_simulation.py FS-cell CONFIG1 -s --SEED 4 --tstop 10 --discret_Fe 10  --discret_Fi 10 --max_Fe 30
# paper's value -> 5-10 minutes sim

from modeling_mesoscopic_dynamics.transfer_functions.theoretical_tools import make_fit_from_data
P_FS = make_fit_from_data('data/FS-cell_CONFIG1.npy', with_square_terms=True, verbose=False)

# now displaying the numerical values and the analytical function
from modeling_mesoscopic_dynamics.transfer_functions.plots import make_exc_inh_fig
make_exc_inh_fig('data/FS-cell_CONFIG1.npy', P=P_FS)

# Note that we also used the network parameters because it requires the synaptic properties (number of synapses)
# they are stored in the following module
from modeling_mesoscopic_dynamics.synapses_and_connectivity.syn_and_connec_library import get_connectivity_and_synapses_matrix
print(get_connectivity_and_synapses_matrix('CONFIG1'))

%run modeling_mesoscopic_dynamics/network_simulations/ntwk_sim_demo.py --CONFIG RS-cell--FS-cell--CONFIG1 -f data/config1.npy

# showing the ouput of the simulations
from modeling_mesoscopic_dynamics.network_simulations.plot_single_sim import plot_ntwk_sim_output
DATA = np.load('data/config1.npy')
AX, FIG = plot_ntwk_sim_output(*DATA,\
                               zoom_conditions=[0,1000], bar_ms=100,\
                               raster_number=500)

# now comparing the simulation with the mean-field
from modeling_mesoscopic_dynamics.network_simulations.compare_with_mean_field import plot_ntwk_sim_output
FIG = plot_ntwk_sim_output(*DATA, min_time=500) # we discard the initial 500ms to evaluate the dynamics
# the fixed point of the mean field is found by launching a trajectory
# of the mean-field dynamical system and taking the point of convergence
# see modeling_mesoscopic_dynamics/mean_field

# running simulation
%run modeling_mesoscopic_dynamics/network_simulations/waveform_input.py --CONFIG RS-cell--FS-cell--CONFIG1 -f data/waveform.npy --sim

# comparing with the time-varying mean-field
# again the dynamical system is implemented in modeling_mesoscopic_dynamics/mean_field
%run modeling_mesoscopic_dynamics/network_simulations/waveform_input.py --CONFIG RS-cell--FS-cell--CONFIG1 -f data/waveform.npy

%run modeling_mesoscopic_dynamics/ring_model/single_trial.py

%run modeling_mesoscopic_dynamics/ring_model/single_trial.py --conduction_velocity_mm_s 50

%run modeling_mesoscopic_dynamics/ring_model/single_trial.py --Tau1 10e-3 --Tau2 50e-3

%run modeling_mesoscopic_dynamics/ring_model/single_trial.py --exc_connect_extent 1 --sX 4

"""
plotting the VSD data and evidence the propagation
"""
from scipy.io import loadmat
import numpy as np
from graphs.my_graph import set_plot
from scipy.signal import convolve2d
import matplotlib.cm as cm

def load_data(filename, Nsmooth=1, tshift=0):
    f = loadmat(filename)
    data = 1e3*f['matNL'][0]['stim1'][0]
    time = f['matNL'][0]['time'][0].flatten()+tshift
    space = f['matNL'][0]['space'][0].flatten()
    if Nsmooth>0:
        smoothing = np.ones((Nsmooth, Nsmooth))/Nsmooth**2
        smooth_data = convolve2d(data, smoothing, mode='same')
    else:
        smooth_data = data
    return time, space, smooth_data
    
def plot_response(filename,
                  tshift=0, Nsmooth=2,
                  t0=-np.inf,t1=np.inf,
                  Nlevels=10,
                  vsd_ticks=[-0.5, 0, 0.5, 1.]):
    
    fig, ax = plt.subplots(1, figsize=(4.7,3))
    plt.subplots_adjust(bottom=.23, top=.9, right=.84, left=.25)

    time, space, smooth_data = load_data(filename, Nsmooth=Nsmooth, tshift=tshift)
    cond = (time>t0) & (time<t1)
    c = ax.contourf(time[cond], space, smooth_data[:,cond],\
             np.linspace(smooth_data.min(), smooth_data.max(), Nlevels), cmap=cm.viridis)
    plt.colorbar(c, label='VSD signal ($\perthousand$)',
                 ticks=vsd_ticks)

    set_plot(ax, xlabel='time (ms)', ylabel='space (mm)')
    return fig, ax

# full response
fig, _ = plot_response('data/VSD_data_session_example.mat')

# zoomed response
fig, ax = plot_response('data/VSD_data_session_example.mat',
                        t0=-50, t1=100, Nsmooth=0)
time, space, data = load_data('data/VSD_data_session_example.mat', Nsmooth=0, tshift=0)

"""
amp_criteria is the threshold of the crossing of maximum amplitude
signal_criteria is a minimum amplitude criteria, beyond this level of the maximum observed signal we consider that there is no real evoked reponse, it does not goes beyond noise level
"""
signal_criteria, amp_criteria = 0.4, 0.6
XX, TT = [], []
for i in range(data.T.shape[1]):
    imax = np.argmax(data.T[:,i])
    if data.T[imax,i]>=signal_criteria:
        ii = np.argmin(np.abs(data.T[:imax,i]-amp_criteria*data.T[imax,i]))
        XX.append(space[i])
        TT.append(time[ii])
# ax.plot(np.array(TT), np.array(XX), 'wo', ms=1)        
tt, xx = np.array(TT), np.array(XX)
# for intervals in [[0,2.3], [2.5,5.7], [5.9,8.5]]:
for intervals in [[3, 5.5], [6,8], [8.2, 12]]:
    cond = (xx>intervals[0]) & (xx<intervals[1]) & (tt<20)
    pol = np.polyfit(xx[cond], tt[cond], 1)
    xxx = np.linspace(xx[cond][0], xx[cond][-1])
    plt.plot(np.polyval(pol, xxx), xxx, 'w-', lw=2)
set_plot(ax, xlabel='time (ms)', ylabel='space (mm)', ylim=[3,11])

# we build up the scan run many configurations over a grid
# for testing
%run modeling_mesoscopic_dynamics/ring_model/parameter_scan.py --N 2 
# UNCOMMENT for paper value
# %run modeling_mesoscopic_dynamics/ring_model/parameter_scan.py --N 5

# this writes a bash file in the 'modeling_mesoscopic_dynamics/ring_model' folder
# launch it with "bash modeling_mesoscopic_dynamics/ring_model/bash_parameter_scan.sh"
# (but ideally on a server, this is a bit long ~5h of simulations)
print('execute the following script')
!(ls modeling_mesoscopic_dynamics/ring_model/bash_parameter_scan.sh)

# we align the data and the model to compute a residual:
from modeling_mesoscopic_dynamics.ring_model.compare_to_model import get_data, get_residual
new_time, space, new_data = get_data('data/VSD_data_session_example.mat', t0 = -50, t1=200)
print('Residual with this configuration = ', get_residual(new_time, space, new_data,
                                                          fn='data/example_data.npy',
                                                          with_plot=True))

# then we compute the residual with respect to that datafile for all the scanned configuration
%run modeling_mesoscopic_dynamics/ring_model/parameter_scan.py --analyze

# and we just take the minimum among all configurations, here the red point
# (note that here it is only for N=2, run N=5 to have it identical to the paper)
%run modeling_mesoscopic_dynamics/ring_model/parameter_scan.py --plot

# this procedure is repeated for all session's data to yield the results of Fig. 9
# (N.B. with the N=5 grid discretization)