生物物理学 2019年度レポート（能瀬担当分）

import numpy as np
import matplotlib.pyplot as plt

import nest

# create Hodgkin-Huxley neuron with delta-shaped synaptic currents.
neuron=nest.Create('hh_psc_alpha')
# create a spike generator
spikegenerator=nest.Create('spike_generator')
# modify spike generation with times 10 and 50 ms
nest.SetStatus(spikegenerator,{'spike_times': [10.,50.]})
# create a voltmeter
voltmeter = nest.Create('voltmeter')
# create add a spikedetector!
spikedetector = nest.Create('spike_detector')
# connect spike generator with a given synaptic specification (weight)
nest.Connect(spikegenerator, neuron, syn_spec={'weight':1e3})
# connect voltmeter to the neuron for measurements
nest.Connect(voltmeter, neuron)
# connect spikedetector to the neuron
nest.Connect(neuron, spikedetector)

### run simulation for 100ms
nest.Simulate(100.)

# read out recording time and voltage from voltmeter (check the parameter list!)
times=nest.GetStatus(voltmeter)[0]['events']['times']
voltage=nest.GetStatus(voltmeter)[0]['events']['V_m']

# plot results
plt.plot(times,voltage)
plt.xlabel('Time (ms)');
plt.ylabel('Membrane potential (mV)')
plt.grid(linestyle='-', linewidth=.25, alpha=.7)
plt.show()

Ek = -80 # mV
gKmax = 83 # Ω^-1/cm^2
Vc = [25, 60, 85, 100] # mV
Ikinf = [25,60,85,100] # mA/cm^2

V = np.linspace(-70,50,1000)
alpha = lambda V: 0.01*(10-V)/(np.exp((10-V)/10)-1)
beta  = lambda V: 0.125*np.exp(-V/80)

plt.title("The relationship between\n'Membrane potential(V)' and 'alpha','beta' (6.3 degrees Celsius)")
plt.plot(V,alpha(V))
plt.plot(V,beta(V))
plt.xlabel("Membrance potential (mV)"), plt.ylabel("alpha,beta")
plt.grid()

n_inf = lambda V:alpha(V)/(alpha(V)+beta(V))

f = np.array([0,0,0,1,1,1])
g = np.array([1,1,0,1,0,0])
A = np.zeros(shape=(g.shape+f.shape))
k = 1

A += k * np.expand_dims(f, axis=0) * np.expand_dims(g, axis=1)

print(A)

G = A.dot(f)

print("G:", G.astype(int))
print("g:", g)

import nest
import time
import random
import numpy as np
import matplotlib.pyplot as plt

def plotPatternRaster(spikes, ax=None, cmap="binary_r", **kwargs):
    if ax is None: fig,ax=plt.subplots()
    H,W = spikes.shape
    ax.matshow(spikes, aspect=W/H, cmap=cmap, **kwargs)
    ax.set_title("Raster plot of presynaptic neurons"), ax.set_xlabel("time (ms)"), ax.set_ylabel("neuron ID")
    return ax

def plotFiringRates(spikes, fig=None):
    if fig is None: fig = plt.figure()
    num_neurons, _ = spikes.shape
    firing_rates = np.mean(spikes, axis=1)*100
    neurons_IDs  = np.arange(num_neurons)
    #=== Plot ===
    ax1 = fig.add_subplot(2,1,1)
    ax1.barh(neurons_IDs, firing_rates, height=1)
    ax1.set_title("Average firing rates of presynaptic neurons"), ax1.set_ylabel("Neuron ID")

    ax2 = fig.add_subplot(2,1,2)
    ax2.hist(firing_rates)
    ax2.set_xlabel("Average firing rate (Hz)"), ax2.set_ylabel("Frequency")
    return fig

class PoissonProcessGenerator():
    def __init__(self, **params):
        if len(params) > 0: self.load_params(**params)
        self.repeated=False
        self.pattern_repeated_spikes = None
        self.simulation_length = None
        self.num_presynaptic_neurons = None

    def load_params(self, **params):
        self.__dict__.update(params)

    def generate_spikes(self, simulation_length=1000, num_presynaptic_neurons=2000, seed=None):
        """Poisson Process.
        - each presynpatic neuron will emit spikes independently
        - generated by a 'Poisson process' which varies randomly in firing rate r in the range of [0,90] Hz
        - the rate of change in r is modified by dr where:
            - started at 0
            - updated randomly from a uniform distribution over [−360,360] Hz/s every millisecond
            - clipped to within [−1800,1800] Hz/s
        - manually add some additional spikes to guarantee that 'in every 50ms time bin, each neuron spikes at least once.'
        ~~~~~~~~~
        Arguments
        @params simulation_length      : (int) the length of the simulation processing. (memorize)
        @params num_presynaptic_neurons: (int) the number of independent presynaptic neurons. (memorize)
        @params seed                   : (int) If you want to make an experiment replicable.
        """
        np.random.seed(seed)
        start = time.time()

        self.simulation_length = simulation_length
        self.num_presynaptic_neurons = num_presynaptic_neurons

        spikes = np.zeros(shape=(num_presynaptic_neurons,simulation_length))
        spikes[:,0] = np.clip(np.random.poisson(self.starting_Hz/1000), 0, 1)
        # ds is randomly picked from a uniform distribution over [−360,+360] Hz/s
        ds = np.random.uniform(self.min_ds, self.max_ds, size=(num_presynaptic_neurons, simulation_length))
        dr = np.zeros_like(ds)
        for t in range(simulation_length-1):
            # dr is clipped to within [−1800,1800] (Hz/s)
            dr[:,t+1] = np.clip(dr[:,t] + ds[:,t], self.min_dr, self.max_dr)
        # r varies beteen [0,90] (Hz)
        r = np.insert(dr[:,:-1], 0, self.starting_Hz, axis=1)
        for t in range(simulation_length-1):
            r[:,t+1] = np.clip(r[:,t] + dr[:,t], self.min_r, self.max_r)
            # Manually add some additional spikes to guarantee that
            # "in every 50ms time bin, each neuron spikes at least once".
            spikes[:,t+1] = np.clip(np.random.poisson(r[:,t+1]/1000), 0, 1)
            if t+1>=self.min_span:
                mask    = np.where(np.sum(spikes[:,t+1-self.min_span:t+1], axis=1)==0)[0]
                add_pos = np.random.randint(t+1-self.min_span, t+1, size=len(mask))
                for idx,pos in zip(mask,add_pos):
                    spikes[idx,pos] = 1
        self.spikes = spikes
        print(f"Processing Time: {time.time()-start:.3f}[s]")
        return spikes

    def select_segment(self, seed=None):
        """Preprocessing.
        1. Divide the presynaptic spiking activity into segments of 50ms (segment_length)
        2. select a random 50ms segment (idx=selection) to be considered as the "pattern"
        """
        self.num_segments = int(self.simulation_length/self.segment_length)
        selection = np.random.RandomState(seed).randint(1, self.num_segments)
        pattern_spikes = self.spikes[:,selection*self.segment_length:(selection+1)*self.segment_length]
        self.pattern_spikes = pattern_spikes
        self.selection = selection
        print(f"Spikes divided into {self.num_segments} segments.")
        print(f"{selection}-th segment was selected as a Pattern.")

    def memorize(self, spikes, repeated=False):
        """ If you already have Poisson Process array (`spont_100s_alt.npy`), you have to memorize it. """
        print("[Before]")
        print(f"simulation length: {self.simulation_length}")
        print(f"Number of the presynaptic neurons: {self.num_presynaptic_neurons}")
        print("Memorizing...")
        self.num_presynaptic_neurons, self.simulation_length = spikes.shape
        if repeated:
            self.pattern_repeated_spikes = spikes
            self.repeated = True
        else:
            self.spikes = spikes
            self.repeated = False
        print("[After]")
        print(f"simulation length: {self.simulation_length}")
        print(f"Number of the presynaptic neurons: {self.num_presynaptic_neurons}")

    def repeat_segment(self, proportion_of_neurons=0.5, frequency_of_patterns=0.25):
        if self.repeated:
            print("Requirement already satisfied.")
            return (None,None)
        else:
            self.pattern_repeated_spikes = np.copy(self.spikes)
            self.repeated_position = np.zeros(self.spikes.shape)
            num_repeats = round(frequency_of_patterns*self.num_segments)
            idx_repeats = random.sample(range(self.num_segments), num_repeats)
            for idx in idx_repeats:
                self.pattern_repeated_spikes[0:int(self.num_presynaptic_neurons*proportion_of_neurons), idx*self.segment_length:(idx+1)*self.segment_length] = self.pattern_spikes[0:int(self.num_presynaptic_neurons*proportion_of_neurons),:]
                self.repeated_position[0:int(self.num_presynaptic_neurons*proportion_of_neurons),idx*self.segment_length:(idx+1)*self.segment_length] = 1
            print(f"Number of spikes added: {np.sum(self.repeated_position)}")
            print(f"Spikes are saved as `self.pattern_repeated_spikes`.")
            self.repeated = True
            return self.pattern_repeated_spikes, self.repeated_position

    def generatePoissonNoise(self, frequency=10):
        noise = np.random.poisson(frequency/1000, size=(self.spikes.shape))
        return noise

    def addPoisonNoise(self, frequency=10):
        spikes = self.pattern_repeated_spikes if self.pattern_repeated_spikes is not None else self.spikes
        noise = self.generatePoissonNoise(frequency=frequency)
        self.noise_added_spikes = np.clip(spikes + noise, 0, 1)
        print(f"Spikes are saved as `self.noise_added_spikes`.")
        return self.noise_added_spikes

class NeuronalNetwork():
    def __init__(self, **params):
        if len(params) > 0: self.load_params(**params)
        self.built = False
    def load_params(self, **params):
        self.__dict__.update(params)
    def build(self, spike_indicies, re=False):
        """build the Neuronal Network
        @param spike_indicies: n-th presynaptic_neuron's i-th spike occurs at spike_indicies[n][i]
        """
        if re: print("Reset Kernel.")
        re = "Re" if re else ""
        nest.ResetKernel()
        self.num_presynaptic_neurons = len(spike_indicies)
        self.spike_generator = nest.Create('spike_generator', self.num_presynaptic_neurons, params=[{"spike_times": np.asarray(spike_idx, dtype=float)+1} for spike_idx in spike_indicies])
        self.parrot_neuron   = nest.Create('parrot_neuron',   self.num_presynaptic_neurons)
        self.post_neuron     = nest.Create('iaf_psc_alpha')
        self.spike_detector  = nest.Create('spike_detector')
        self.voltmeter       = nest.Create('voltmeter', params={'interval': 0.1})
        self.spike_indicies  = spike_indicies
        print(f"{re}Created {self.num_presynaptic_neurons} spike generators.")
        print(f"{re}Created {self.num_presynaptic_neurons} parrot neurons.")
        print(f"{re}Created 1 post_neuron.")
        print(f"{re}Created spike detector.")
        print(f"{re}Created voltmeter.")

    def addGaussianPackets(self, simulation_length_pulses,
                                 proportion_of_neurons_pulsing,
                                 number_of_pulses,
                                 spikes_per_pulse,
                                 std_spikes_per_pulse,
                                 seed=None):
        pulse_center_times = np.random.RandomState(seed).randint(10, simulation_length_pulses-10, number_of_pulses).astype(float)
        self.gaussian_pulses = nest.Create("pulsepacket_generator", params={
            "pulse_times": pulse_center_times,
            "activity": spikes_per_pulse,
            "sdev": std_spikes_per_pulse
        })
        conn_dict = {
            'rule': 'fixed_total_number',
            'N': int(self.num_presynaptic_neurons/proportion_of_neurons_pulsing)
        }
        nest.Connect(self.gaussian_pulses, self.parrot_neuron, conn_dict)

    def compile(self, synapse_model="static", ex=1, **params):
        if self.built: self.build(self.spike_indicies, re=True)
        if len(params) > 0: self.load_params(**params)
        nest.SetStatus(self.post_neuron, [{"tau_minus": self.tau_LTD}])
        nest.Connect(self.spike_generator, self.parrot_neuron, "one_to_one")
        syn_dict = self.SynapseModelHandler(name=synapse_model, ex=ex)
        nest.Connect(self.parrot_neuron, self.post_neuron, "all_to_all", syn_dict)
        nest.Connect(self.post_neuron, self.spike_detector)
        nest.Connect(self.voltmeter, self.post_neuron)
        print(f"Connected {self.num_presynaptic_neurons} spike generators to {self.num_presynaptic_neurons} parrot neurons. (one-to-one)")
        print(f"Connected {self.num_presynaptic_neurons} parrot neurons to 1 post neuron. (all-to-all) Method: {syn_dict['model']}")
        print("Connected post neuron to spike detector")
        print("Connected voltmeter to post neuron")
        self.built = True

    def simulate(self, T=1000., ax=None, **plotkwargs):
        if ax is None: fig, ax = plt.subplots()
        nest.Simulate(float(T))

        times      = nest.GetStatus(self.voltmeter)[0]['events']['times']
        voltage    = nest.GetStatus(self.voltmeter)[0]['events']['V_m']
        spiketimes = nest.GetStatus(self.spike_detector)[0]['events']['times']

        # As the simulation span is larger than reaction length,
        # we will modify the voltage to Make it understandable.
        for spike in spiketimes:
            voltage[(np.abs(np.asarray(times) - spike)).argmin()] = 0.

        ax.plot(times,voltage, **plotkwargs)
        ax.set_xlabel('Time (ms)'), ax.set_ylabel('Membrane potential (mV)'), ax.grid(linestyle='-', linewidth=.25, alpha=.7)
        return ax

    def mask_pattern(self, ax, is_pattern, simulation_times):
        add_minus_1 = np.insert(arr=np.any(is_pattern[:,:simulation_times], axis=0), obj=0, values=False).astype(int)
        masks = add_minus_1[1:] - add_minus_1[:-1]
        start_points = np.where(masks == 1)[0]
        end_points = np.where(masks == -1)[0]
        for s,e in zip(start_points, end_points):
            ax.axvspan(s, e, color = "coral", alpha=0.7)
        return ax

    def plotWeights(self, is_pattern_neuron, ax=None):
        if ax is None:
            fig, ax = plt.subplots()
        pre_to_post_conninfo = nest.GetConnections(self.parrot_neuron, self.post_neuron)
        pre_to_post_weights  = np.asarray(nest.GetStatus(pre_to_post_conninfo, keys='weight'))
        idxes = np.arange(len(pre_to_post_weights))
        ax.barh(idxes[~is_pattern_neuron], pre_to_post_weights[~is_pattern_neuron], color="blue", label="Random Neuron")
        ax.barh(idxes[is_pattern_neuron], pre_to_post_weights[is_pattern_neuron],    color="red", label="Pattern Neuron")
        ax.set_xlabel("Synaptic weights"), ax.set_ylabel("Neuron ID."), ax.legend()
        return ax

    def SynapseModelHandler(self,name,ex):
        if ex==1:
            SynapseModelDict = {
                "static"  : {
                    "model": "static_synapse",
                    "weight": 0.5,
                },
                "standard": {
                    "model": "stdp_synapse",
                    "weight": 1.0,
                    "tau_plus": self.tau_LTP,
                    "mu_plus": self.a_LTP,
                    "mu_minus": self.a_LTD,
                },
                "triplet" : {
                    "model": "stdp_triplet_synapse",
                    "weight": 1.0,
                    "tau_plus": self.tau_LTP,
                    "Aminus": self.a_LTD/2,
                    "Aminus_triplet": self.a_LTD/10,
                    "Aplus": self.a_LTP,
                    "Aplus_triplet": self.a_LTP*10,
                 },
                 "FACETS": {
                    "model": "stdp_facetshw_synapse_hom",
                    "weight": 0.7,
                    "tau_plus": self.tau_LTP,
                    "tau_minus_stdp": self.tau_LTD,
                 }
            }
            if name not in SynapseModelDict:
                raise ValueError(f"Couldn't understand {name}. Please chose from\033[01m\033[36m '{', '.join(SynapseModelDict.keys())}' \033[0m")
            return SynapseModelDict[name]
        elif ex==2:
            SynapseModelDict = {
                "static"  : {
                    "model": "static_synapse",
                    "weight": 1.0,
                },
                "standard": {
                    "model": "stdp_synapse",
                    "weight": 1.0,
                    "tau_plus": self.tau_LTP,
                    "mu_plus": self.a_LTP*1000,
                    "mu_minus": self.a_LTD/1000,
                    "Wmax": 5.0
                },
                "triplet" : {
                    "model": "stdp_triplet_synapse",
                    "weight": 1.0,
                    "tau_plus": self.tau_LTP,
                    "Aminus": self.a_LTD/2,
                    "Aminus_triplet": self.a_LTD/10,
                    "Aplus": self.a_LTP,
                    "Aplus_triplet": self.a_LTP*10,
                    "Wmax": 5.0
                 },
                 "FACETS": {
                    "model": "stdp_facetshw_synapse_hom",
                    "weight": 0.75,
                    "tau_plus": self.tau_LTP,
                    "tau_minus_stdp": self.tau_LTD,
                    "Wmax": 5.0
                 }
            }
            if name not in SynapseModelDict:
                raise ValueError(f"Couldn't understand {name}. Please chose from\033[01m\033[36m '{', '.join(SynapseModelDict.keys())}' \033[0m")
            return SynapseModelDict[name]

poisson_params = {
    "min_span": 50, # [ms]
    
    "starting_Hz": 54, # [Hz]
    "min_r": 0,        # [Hz]
    "max_r": 90,       # [Hz]
    
    "min_ds": -3.6e-1, # [Hz/ms]
    "max_ds": 3.6e-1,  # [Hz/ms]
    "min_dr": -1.8,    # [Hz/ms]
    "max_dr": 1.8,     # [Hz/ms]
    
    "segment_length": 50, # [ms] Length of 1 Presynaptic spiking segment length.
}

seed = 0

num_presynaptic_neurons = 2000 # Number of the presynaptic neurons.
simulation_length = 1000       # [ms]

poisson_generator = PoissonProcessGenerator()
poisson_generator.load_params(**poisson_params)

spikes = poisson_generator.generate_spikes(simulation_length, num_presynaptic_neurons, seed=seed)

poisson_generator.select_segment(seed=seed)

fig,ax = plt.subplots(figsize=(5,5))
ax = plotPatternRaster(spikes, ax=ax)
plt.tight_layout()
plt.show()

proportion_of_neurons=0.5
frequency_of_patterns=0.25

pre_spikes_wrepeats, output_identities = poisson_generator.repeat_segment(proportion_of_neurons, frequency_of_patterns)

frequency = 10 # [Hz]

noise_added_spikes = poisson_generator.addPoisonNoise(frequency=frequency)

fig,ax = plt.subplots(figsize=(5,5))
ax = plotPatternRaster(noise_added_spikes, ax=ax)
plt.tight_layout()
plt.show()

fig,ax = plt.subplots(figsize=(5,5))
ax = plotPatternRaster(noise_added_spikes, ax=ax)
ax = plotPatternRaster(output_identities, cmap="RdGy_r", ax=ax, alpha=0.7)
plt.tight_layout()
plt.show()

simulation_length_100s = 100_000
num_presynaptic_neurons = 2_000
poisson_generator = PoissonProcessGenerator()
poisson_generator.load_params(**poisson_params)
pre_spikes_wrepeats_100s = poisson_generator.generate_spikes(simulation_length_100s, num_presynaptic_neurons, seed=seed)

proportion_of_neurons = 0.5
frequency_of_patterns = 0.25

poisson_generator.select_segment(seed=seed)
pre_spikes_wrepeats, output_identities = poisson_generator.repeat_segment(proportion_of_neurons, frequency_of_patterns)
noise_added_spikes = poisson_generator.addPoisonNoise(frequency=frequency)

repeated_idxes = np.any(poisson_generator.repeated_position, axis=1)

print(f"Pattern shape = {poisson_generator.pattern_spikes.shape}")
print(f"Poisson Process shape = {noise_added_spikes.shape}")

simulation_times = 1000

spike_indicies = [np.where(is_spike)[0] for is_spike in noise_added_spikes==1]

STDP_paramaters = {
    'a_LTP' : 0.03125,
    'a_LTD' : 0.85 * 0.03125,
    'tau_LTP' : 16.8,
    'tau_LTD' : 33.7,
}

model = NeuronalNetwork(**STDP_paramaters)
model.build(spike_indicies)

model.compile(synapse_model="static", ex=1)

ax = model.simulate(T=simulation_times, color="black")
ax = model.mask_pattern(ax, output_identities, simulation_times)
plt.show()

ax = model.plotWeights(repeated_idxes)
plt.show()

model.compile(synapse_model="standard", ex=1)

ax = model.simulate(T=simulation_times, color="black")
ax = model.mask_pattern(ax, output_identities, simulation_times)
plt.show()

ax = model.plotWeights(repeated_idxes)
plt.show()

model.compile(synapse_model="triplet", ex=1)

ax = model.simulate(T=simulation_times, color="black")
ax = model.mask_pattern(ax, output_identities, simulation_times)
plt.show()

ax = model.plotWeights(repeated_idxes)
plt.show()

model.compile(synapse_model="FACETS", ex=1)

ax = model.simulate(T=simulation_times, color="black")
ax = model.mask_pattern(ax, output_identities, simulation_times)
plt.show()

ax = model.plotWeights(repeated_idxes)
plt.show()