使用 solve_ivp 进行同时事件检测

Question

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

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

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

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

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

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

这也捕获了图中的双重事件和结果

使用 solve_ivp 进行同时事件检测

Simultaneous event detection using solve_ivp

python

scipy

ode