Source code for emodelrunner.synplas_analysis

"""Analysis tools for synapse plasticity."""

# Copyright 2020-2022 Blue Brain Project / EPFL

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at


# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

import h5py
import numpy as np

[docs] def get_epsp_vector(t, v, spikes, window): """Extract EPSPs at time `spikes` from voltage trace `v`. Args: t (numpy.ndarray): Time vector (ms). v (numpy.ndarray): Soma voltage trace. Must have the same size of `t` (mV). spikes (numpy.ndarray): Time of presynaptic spikes. Every value in `spikes` must be in the interval [min(`t`), max(`t`)] (ms). window (float): Size of EPS detection window in ms. Returns: numpy.ndarray: Vector of EPSPs (mV). Raises: RuntimeError: if a postsynaptic cell spike is found during a connectivity test RuntimeError: if the detection window can have multiple EPSPs in it """ # Verify presence of only one EPSP in window if np.max(np.diff(spikes)) < window: raise RuntimeError("Only one EPSP should be present in the detection window") # Get EPSPs n = len(spikes) epsps = np.zeros(n) for i in range(n): w0 = np.searchsorted(t, spikes[i]) # Beginning of EPSP window vbaseline = v[w0] w1 = np.searchsorted(t, spikes[i] + window) # End of EPSP window epspmaxidx = np.argmax(v[w0:w1]) + w0 vepsp = v[epspmaxidx] # Abort if postsynaptic spike if vepsp > -30: raise RuntimeError("Postsynaptic cell spiking during connectivity test") # Compute epsps epsp = vepsp - vbaseline epsps[i] = epsp return epsps
[docs] def epsp_slope(vtrace): """Returns slope between points at [0.3, 0.75] * peak of given trace. Attention! This function expects that the trace is sampled (or interpolated) at a single timestep, since time is not taken into account during calculations. Args: vtrace (numpy.ndarray): interpolated voltage trace of one EPSP (mV) Returns: float: slope of the EPSP (mV) """ # Remove baseline v = vtrace - vtrace[0] # Find peak peak = np.max(v) peak_idx = np.argmax(v) # Get 30% and 75% rise time indices idx0 = np.argmin(np.abs(v[:peak_idx] - 0.3 * peak)) idx1 = np.argmin(np.abs(v[:peak_idx] - 0.75 * peak)) # Get slope of fitting line between the two m = (v[idx1] - v[idx0]) / (idx1 - idx0) return m
[docs] class Experiment: """A full STDP induction experiment. The experiment consists of two connectivity tests (C01 and C02), separed by the induction protocol. Attributes: t (numpy.ndarray): Time vector (ms) v (numpy.ndarray): Soma voltage trace (ms) spikes (numpy.ndarray): Pre-synaptic spike times (ms) duration (dict): Duration of C01 and C02 phases (ms) period (float): The Cx period (ms) c01period (float): The C01 period (ms) c02period (float): The C02 period (ms) epspwindow (float): Size of EPS detection window (ms) cxs (list): Names of the Cx tests. These names are used as keys in the dictionaries (cxspikes, duration, epsp, cxtrace) cxspikes (dict): Pre-synaptic spikes times (ms) occuring during each test (C01 & C02) """ def __init__( self, data, c01duration=40.0, c02duration=40.0, period=10.0, c01period=None, c02period=None, ): """Constructor. Args: data (str or dict): Path to simulation results file or results dictionary c01duration (float): Duration of C01 phase in minutes. c02duration (float): Duration of C02 phase in minutes. period (float): The Cx period in seconds. c01period (float): The C01 period in seconds. If None, period will be used. c02period (float): The C02 period in seconds. If None, period will be used. Raises: TypeError if data is neither a dictionary nor a string """ if isinstance(data, dict): self.t = data["t"] self.v = data["v"] self.spikes = data["prespikes"] elif isinstance(data, str): h5file = h5py.File(data, "r") self.t = h5file["t"][()] self.v = h5file["v"][()] self.spikes = h5file["prespikes"][()] h5file.close() else: raise TypeError("data argument of Experiment has to be dict or str") # Store other attributes self.duration = { "C01": c01duration * 60 * 1000.0, # min to ms "C02": c02duration * 60 * 1000.0, } # min to ms self.period = period * 1000.0 # sec to ms self.c01period = ( c01period * 1000.0 if c01period is not None else period * 1000.0 ) self.c02period = ( c02period * 1000.0 if c02period is not None else period * 1000.0 ) # Create common attributes self.epspwindow = 100.0 # ms self.cxs = ["C01", "C02"] self.cxspikes = { "C01": self.spikes[: int(self.duration["C01"] / self.c01period)], "C02": self.spikes[-int(self.duration["C02"] / self.c02period) :], } @property def epsp(self): """Returns dict containing EPSP vectors for each test (mV).""" epsp = { cx: get_epsp_vector(self.t, self.v, self.cxspikes[cx], self.epspwindow) for cx in self.cxs } return epsp @property def cxtrace(self): """Returns time and interpolated trace of each EPSP as dict.""" tdense = np.linspace(0, self.epspwindow, int(self.epspwindow / 0.025)) cxtrace = {"t": tdense} for cx in self.cxs: cxtrace[cx] = [] for s in self.cxspikes[cx]: idx0 = np.searchsorted(self.t, s) idx1 = np.searchsorted(self.t, s + self.epspwindow) vdense = np.interp(tdense + s, self.t[idx0:idx1], self.v[idx0:idx1]) cxtrace[cx].append(vdense) cxtrace[cx] = np.array(cxtrace[cx]) return cxtrace
[docs] def compute_epsp_interval(self, interval): """Compute mean EPSP amplitude at regular intervals. Args: interval (float): The interval in minutes. Returns: dict: A dictionary of interval statistics for each Cx in Experiment. """ results = {} for cx in self.cxs: n = int(self.duration[cx] / (interval * 60 * 1000.0)) epsp_groups = np.split(self.epsp[cx], n) spike_groups = np.split(self.cxspikes[cx], n) avg = np.mean(epsp_groups, axis=1) sem = np.std(epsp_groups, axis=1) / np.sqrt(len(epsp_groups[0])) t = np.mean(spike_groups, axis=1) results[cx] = {"avg": avg, "sem": sem, "t": t} return results
[docs] def compute_epsp_ratio(self, n, method="amplitude", full=False): """Compute mean EPSP change. Args: n (int): Number of sweeps in Cx to be considered for mean EPSP calculation method (str): Method used to compute EPSP ratio (amplitude or slope) full (bool): whether to return the mean value and std of epsp before and after cell stimulus in addition to the epsp ratio Returns: if full is False, returns - ratio_epsp (float): Mean EPSP change if full is True, returns a tuple containing - epsp_before (float): mean value of EPSP or slope of mean EPSP trace in test C01 - epsp_after (float): mean value of EPSP or slope of mean EPSP trace in test C02 - ratio_epsp (float): Mean EPSP change - epsp_before_std (float): std of EPSP in test C01 if method is amplitude, else 0 - epsp_after_std (float): std of EPSP in test C02 if method is amplitude, else 0 Raises: ValueError: if method is neither 'amplitude' nor 'slope' """ if method == "amplitude": epsp_before = np.mean(self.epsp["C01"][-n:]) epsp_before_std = np.std(self.epsp["C01"][-n:]) epsp_after = np.mean(self.epsp["C02"][-n:]) epsp_after_std = np.std(self.epsp["C02"][-n:]) elif method == "slope": epsp_before = epsp_slope(np.mean(self.cxtrace["C01"][-n:], axis=0)) epsp_before_std = 0 epsp_after = epsp_slope(np.mean(self.cxtrace["C02"][-n:], axis=0)) epsp_after_std = 0 else: raise ValueError(f"Unknown method {method}") epsp_ratio = epsp_after / epsp_before if full: return epsp_before, epsp_after, epsp_ratio, epsp_before_std, epsp_after_std else: return epsp_ratio
[docs] def normalize_time(self, t): """Normalize time vector `t`. Convert milliseconds to minutes and shift t = 0 to beginning of induction phase. Args: t (numpy.ndarray): Time vector (ms) Results: tnorm (numpy.ndarray): Normalized time vector (mn) """ return (t - self.duration["C01"]) / (60.0 * 1000.0)