[FIXED] How to react to simultaneous threshold crossings in solve_ivp in Python

Issue

I am considering a neural system in which several segments are connected in one dimension and the time derivative of the activity of a segment is affected by the activity of the neighboring segments. In other words, let the activity of the i-th segment be E_i, which can be expressed as

dE_i/dt = f(a E_i(t) + b_{i-1}(t)E_{i-1}(t) + c_i(t)E_{i+1}(t)))

using a function f. The b and c represent the synaptic strength between segments. a is a constant, but b and c are time-varying. If the activity of the two segments across the synapse exceeds the threshold at the same time, b and c will be modified.

To simulate such a model, I coded the following.

import numpy as np
from scipy.integrate import solve_ivp

seg = 8

t0 =
tend = 
dt = 
threshold = 

x0 = [0]*seg

def event(t, x):
　　?

def model(t, x):
　　・
　　・
　　・
　　return (dE)

While True:
　　t_step = np.arange(t0, tend - dt, dt)
　　solve = solve_ivp(model, [t0, tend], x0, t_eval = t_step, events = event)

　　if solve.status == 1:
　　　・
　　　・
　　　・

It solves a simultaneous differential equation using solve_ivp, and when there is an i (i=0 to seg-1) where E_i and E_{i+1} exceed the threshold at the same time, the solver pauses, modifies b and c appropriately, and solves the differential equation again. I've just started learning Python, so I don't know anything about it, but I would appreciate it if you could tell me how to do it.

※ "Exceeding the threshold at the same time" means that while one segment is higher than the threshold, the neighboring segment exceeds the threshold. It does not mean that two segments cross the threshold at exactly the same time.
I'm sorry for my poor code explanation and my poor English.

Solution

You get close to the limitations of the event mechanism here. What you have to do is implement one event function per synapse, and pass the array of events to solve_ivp. This can be done via lambda functions

def event_proto(i,t,E): return min(E[i],E[i+1])-threshold

events = [ lambda t,E: event_proto(i,t,E) for i in range(seg-1) ]

I'm not sure if you have to set the direction and terminal condition once on the prototype function or separately on all functions of the array. You can find out which event triggered from solve.t_events. Also check if several events can trigger on the same time segment, then you would have to determine the minimal time for the first event.

Answered By - Lutz Lehmann