使用solve_ivp进行同时事件检测

Question

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

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

我做了一些基本的分析，看解数组，我的猜测是scipy.solve_ivp的事件检测无法检测到两个神经元何时触发。我这么说的原因是，第 2 和第 3 神经元的放电率应该相同，因为它们只有具有外部电流的神经元是第 1 神经元。

但是，当它们同时触发时，事件检测仅检测到一个事件，因此无法重置另一个事件的电压，因此呈指数增长。

我的最终目标是将其与其他类型的神经元耦合，但由于其中许多具有内在的复极化动力学，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，但由于 deadlines，我想让这个脚本正常工作，因为即使是 QIF 神经网络的快速而肮脏的实现网络会做我的初步分析。

我是一名生物学专业的学生，​​只知道一点 Python 和 MATLAB，但无论如何我都将不胜感激。

Answer 1

您确实是正确的，solve_ivp 没有检测到同时发生的其他事件（在您复制组件的情况之外，因为在这里它不太可能在数值模拟中达到这种情况）。您可以手动测试它，因为事件是事件函数的根。所以设置

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]}')

这也会在图中捕获双重事件和结果