Одновременное обнаружение событий с использованием solve_ivp - PullRequest
Новая
       14

Одновременное обнаружение событий с использованием solve_ivp

2 голосов
/ 09 февраля 2020

Я пытаюсь соединить несколько Quadrati c встроенных-и-огненных нейронов.

Мой скрипт успешно работает с двумя нейронами, но когда я изменил скрипт для 3 нейронов Я заметил, что напряжение третьего нейрона внезапно взрывается и, следовательно, интеграция не удалась.

Я провел некоторый базовый c анализ, и, глядя на массив решений, я предполагаю, что обнаружение событий scipy.solve_ivp не может определить, когда два нейрона срабатывают одновременно. Моя причина в том, что частота срабатывания 2-го и 3-го нейронов должна быть одинаковой, поскольку 1-й из них - это только нейрон с внешним током.

Однако, хотя они оба срабатывают вместе, обнаружение только событий обнаруживает одно событие и, следовательно, не может сбросить напряжение другого, отсюда и экспоненциальный рост.

Моя конечная цель состояла бы в том, чтобы связать это с другими типами нейронов, но так как многие из них имеют встроенные c динамика реполяризации, обработка событий QIF - важнейшая часть масштабирования сети. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Define vectors, indices and parameters
resetV = -0.1
nIN = 3
incIN = nIN
ylen = nIN*(incIN)
indIN = np.arange(0,ylen,incIN)
INs = np.arange(0,nIN)    
gI = -0.4    
Ileak = 0.5


# Define heaviside function for synaptic gates (just a continuous step function)
def heaviside(v,thresh):
    H =  0.5*(1 +np.tanh((v-thresh)/1e-8))
    return H

# Define event functions and set them as terminal
def event(t, y):
    return y[indIN[0]] - 2
event.terminal = True

def event2(t,y):
    return y[indIN[1]] - 2
event2.terminal = True

def event3(t,y):
    return y[indIN[2]] - 2
event3.terminal = True

#ODE function
def Network(t,y):

    V1 = y[0]
    n11 = y[1]
    n12 = y[2]
    V2 = y[3]
    n21 = y[4]
    n22 = y[5]
    V3 = y[6]
    n31 = y[7]
    n32 = y[8]

    H = heaviside(np.array([V1,V2,V3]),INthresh)

    dydt = [V1*V1 - gI*n11*(V2)- gI*n12*(V3)+0.5,
            H[1]*5*(1-n11) - (0.9*n11),
            H[2]*5*(1-n12) - (0.9*n12),
            V2*V2 -gI*n21*(V1)- gI*n22*(V3),
            H[0]*5*(1-n21) - (0.9*n21),
            H[2]*5*(1-n22) - (0.9*n22),
            V3*V3 -gI*n31*(V1)- gI*n32*(V2),
            H[0]*5*(1-n31) - (0.9*n31),
            H[1]*5*(1-n32) - (0.9*n32)
            ]

    return dydt

# Preallocation of some vectors (mostly not crucial)
INthresh = 0.5
dydt = [0]*ylen
INheavies = np.zeros((nIN,))
preInhVs = np.zeros((nIN,))          
y = np.zeros((ylen,))

allt = []
ally = []
t = 0
end = 100

# Integrate until an event is hit, reset the spikes, and use the last time step and y-value to continue integration

while True:

    net = solve_ivp(Network, (t, end), y, events= [event,event2,event3])
    allt.append(net.t)
    ally.append(net.y)
    if net.status == 1:

        t = net.t[-1]
        y = net.y[:, -1].copy()
        for i in INs:    
            if net.t_events[i].size != 0:
                y[indIN[i]] = resetV
                print('reseting V%d' %(i+1))

    elif net.status == -1:
        print('failed!')
        print(y[0])
        break
    else:
        break

# Putting things together and plotting
Tp = np.concatenate(ts)
Yp = np.concatenate(ys, axis=1)
fig = plt.figure(facecolor='w', edgecolor='k')
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)
ax1.plot(Tp, Yp[0].T)
ax2.plot(Tp, Yp[3].T)
ax3.plot(Tp, Yp[6].T)
plt.subplots_adjust(hspace=0.8)

plt.show()

Конечно, это только предположение.

Я сейчас нахожусь желая научиться работать с PyDSTool , но из-за крайних сроков я бы хотел, чтобы этот скрипт работал, поскольку даже быстрая и грязная реализация нейронной сети QIF могла бы мой предварительный анализ.

Я изучаю биологию и знаю немного Python и MATLAB, но я был бы признателен за любой вклад, независимо от того. * 10 28 *

1 Ответ

1 голос
/ 09 февраля 2020

Вы действительно правы, solve_ivp не обнаруживает дополнительные события, которые происходят в одно и то же время (за исключением ситуаций, когда вы дублируете компонент, поскольку здесь очень маловероятно, что такая ситуация возникнет при численном моделировании). Вы можете проверить это вручную, поскольку событие является root функции события. Поэтому установите 

def gen_event(i):    
    def event(t, y):
        return y[indIN[i]] - 2
    event.terminal = True
    return event

events = [gen_event(i) for i in range(3)]

и замените тест, для которого функция вызывает событие, на 

    t = net.t[-1]
    y = net.y[:, -1].copy()
    for i in INs:    
        if abs(events[i](t,y)) < 1e-12:
            y[indIN[i]] = resetV
            print(f'reseting V{i+1} at time {net.t_events[i]}')

. Затем он также захватывает двойные события и приводит к графикам

pulses and helper functions

...