Spike-Triggered Average：在神经元放电之前

import numpy as npimport matplotlib.pyplot as plt
# Seed for reproducibilitynp.random.seed(42)
# Parameters for synthetic datasampling_rate = 1000  # in Hzduration = 2  # in secondstime = np.arange(0, duration, 1/sampling_rate)  # time vector

frequencies = [3, 7, 13]  # multiple frequencies for more complexitystimulus = sum(np.sin(2 * np.pi * f * time) for f in frequencies)  # sum of sinusoidsnoise = 0.3 * np.random.randn(len(time))  # increased noise levelstimulus += noisestimulus *= 0.5  # scale down to match example range

# Plot stimulus potential immediately after generating the dataplt.figure(figsize=(6, 3))plt.plot(time, stimulus, 'k')plt.title('Stimulus Potential')plt.ylabel('Potential (mV)')plt.xlabel('Time (s)')plt.ylim(-2, 2)plt.grid(True)plt.show()

threshold = 0.6  # threshold for spike generation, modify for realistic firingspike_probability = np.clip((stimulus - threshold), 0, None)  # only positive values above thresholdspike_probability /= np.max(spike_probability)  # normalize to max probability 1spike_probability *= 0.4  # further scale down to make spikes even sparser
# Convert spike probability to a rate (spikes per second)spike_rate = spike_probability * sampling_rate  # rate in spikes per second
# Generate spike train using Poisson processspikes = np.zeros(len(time))for i in range(len(time)):    spikes[i] = np.random.poisson(spike_rate[i] / sampling_rate)  # Poisson process with rate

# Plot spike train immediately after generating the dataplt.figure(figsize=(6, 3))plt.vlines(time[spikes > 0], ymin=0, ymax=1, color='k')  # Draw vertical lines for spikesplt.title('Evoked Spike Train')plt.xlabel('Time (s)')plt.ylim(-0.1, 1.1)  # Set limits to ensure consistent appearance of spikesplt.yticks([])  # Remove y-axis ticksplt.xticks(np.arange(0, duration + 0.5, 0.5))  # Show x-axis ticks at intervals of 0.5 secondsplt.axhline(y=0, color='k', linestyle='--')  # Add horizontal line at y=0plt.grid(True, which='both', linestyle='--', linewidth=0.5)  # Add gridlines for clarityplt.show()

# Calculate spike-triggered average (STA)window = int(0.3 * sampling_rate)  # 300 ms window for STAsta = np.zeros(window)  # initialize STA for negative time only
# Find indices of spikesspike_indices = np.where(spikes == 1)[0]for idx in spike_indices:    if idx > window:  # only consider past data, not future        sta += stimulus[idx-window:idx]
# Normalize STAif len(spike_indices) > 0:    sta /= len(spike_indices)sta_time = np.arange(-window, 0) / sampling_rate * 1000  # in ms, negative time only

# Plot STA immediately after generating the dataplt.figure(figsize=(6, 4))plt.plot(sta_time, sta * 1000, 'k')plt.title('Spike-Triggered Average Stimulus')plt.xlabel('Time (ms)')plt.ylabel('C (mV)')plt.grid(True)plt.xlim(-300, 0)  # negative time onlyplt.show()

# Generate a new fluctuating stimulus potential with noisefrequencies2 = [2, 7, 10]  # multiple frequencies for more complexitystimulus2 = sum(np.sin(2 * np.pi * f * time) for f in frequencies)  # sum of sinusoidsnoise2 = 0.2 * np.random.randn(len(time))  # increased noise levelstimulus2 += noisestimulus2 *= 0.5  # scale down to match example range

from scipy.stats import zscoredef compute_and_plot_firing_rate(stimulus, sta, time, mode, title):    # Perform convolution to estimate firing rate    predicted_firing_rate = np.convolve(stimulus, sta, mode)
    # Normalize using Z-score    predicted_firing_rate = zscore(predicted_firing_rate)
    # Ensure firing rate is positive    predicted_firing_rate = np.clip(predicted_firing_rate, 0, None)
    # Adjust time vector if mode is 'valid'    if mode == 'valid':        valid_length = len(predicted_firing_rate)        time = time[(len(time) - valid_length) // 2: (len(time) + valid_length) // 2]
    # Plot normalized predicted firing rate    plt.figure(figsize=(6, 4))    plt.plot(time, predicted_firing_rate, 'k')    plt.title(title)    plt.xlabel('Time (s)')    plt.ylabel('Firing Rate')    plt.grid(True)    plt.show()
compute_and_plot_firing_rate(stimulus2, sta, time, 'same', 'Predicted Firing Rate (Same Mode)')compute_and_plot_firing_rate(stimulus, sta, time, 'valid', 'Predicted Firing Rate (Valid Mode)')