使用solve_ivp的同时事件检测>> [

Question

我正在尝试耦合几个Quadratic integrate-and-fire神经元。

我的脚本成功地使用了两个神经元，但是当我修改3个神经元的脚本时，我注意到第三个神经元的电压突然爆炸，因此集成失败。

我做了一些基本分析，然后看解数组，我的猜测是scipy.solve_ivp的事件检测无法同时检测到两个神经元

fire

的时间。我之所以这样说是因为第二和第三神经元的放电速率应该相同，因为它们只有具有外部电流的神经元才是第一电流。但是，虽然它们都一起触发，但事件检测仅检测到一个事件，因此无法重置另一个事件的电压，因此呈指数增长。

我的最终目标是将其与其他类型的神经元耦合，但是由于许多神经元具有内在的复极化动力学，因此对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，但是无论如何我都会很感激。

我正在尝试耦合几个二次积分并发射神经元。我的脚本成功地使用了两个神经元，但是当我修改3个神经元的脚本时，我注意到第三个...

Answer 1

[您的确是正确的，solve_ivp不会检测到同时发生的其他事件（在您复制组件的情况之外，因为在这里很难在数值模拟中达到这种情况）。您可以手动测试，因为事件是事件函数的根。如此设置