基于solve_-ivp的同时事件检测

Question

我试着把几个神经元连接起来

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

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

然而，当它们一起触发时，事件检测只检测到一个事件，因此无法重置另一个事件的电压，从而导致指数增长

我的最终目标是将其与其他类型的神经元结合，但由于其中许多神经元具有内在的复极动力学，QIFs的事件处理是扩展网络的关键部分

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

然后，它还将捕获双事件并在图中显示结果