Блог

Я работаю над моделированием нейронных цепей. Схема состоит из восьми сегментов, и я решаю активность каждого сегмента с помощью дифференциальных уравнений. Дифференциальное уравнение является

dE[i]/dt = f(b[i] * E[i]) (i = 0, 1, ... 7)

и b[i] изменяется при E[i] снижении порога снизу вверх. Чтобы реализовать это, я создал следующий код, используя scipy’s solve_ivp .

 seg_number = 8
t0 = 0
tend = 2e5
x0 = [0]*seg_number

def gen_event(i):    
    def event(t, x):
        return x[i] - threshold
    event.terminal = True
    event.direction = 1
    return event

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


def func(t, x): #differential equation
    dx = []
    for i in range(seg_number):
    ・・・
    return(dx)


while True:
    orbit = solve_ivp(func, (t0, tend), x0, events = events)
    if orbit.status == 1:
        ・・・  # update b[i]

        t0 = orbit.t[-1]
        x0 = orbit.y[:, -1].copy()

    else:
        break
 

Однако в этом коде при b[i] изменении начальных условий и повторном решении дифференциальных уравнений начальное условие определяется как событие, вызывающее бесконечный цикл. Если я добавлю небольшую сумму к начальным условиям, решение продолжится, но я не думаю, что это хорошее решение.

 x0 = orbit.y[:, -1].copy()   1e-12