I would be really thankful for some suggestions how to introduce a delay in one of my ODE in the system and solve it with solve_ivp scipy.

So, I have a system of equations with some analytical functions f1(t), f2(t), f3(t), f4(t) that I pass to solve_ivp function.

def model_step(t,
           y,
           w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9, w_10,
           tau_0, tau_1, tau_2, tau_3,
           threshold,
           power,
           q,
           i_0, i_1, i_2, i_3
           ):
"""
Builds differential equation model
"""

        
f_e, f_p, f_s, f_v = y

ff_e = (min(max((i_0 + w_0 * f_e + w_1 * f1(t) + w_2 * f_p + w_3 * f_s), threshold), 25))
ff_p = (min(max((i_1 + w_4 * f_e + f2(t)), threshold), 25))
ff_s = (min(max((i_2 + w_5 * f_e + w_6 * f3(t) - w_7 * f_v), threshold), 25))
ff_v = (min(max((i_3 + w_8 * f_e - w_9 * f_s +  w_10 * f4(t)), threshold), 25))

    
df_edt = ((q * ff_e ** power) - f_e) / tau_0
df_pdt = ((q * ff_p ** power) - f_p) / tau_1
df_sdt = ((q * ff_s ** power) - f_s) / tau_2
df_vdt = ((q * ff_v ** power) - f_v) / tau_3


dydt = [df_edt, df_pdt, df_sdt, df_vdt]

return dydt

I would like to add a delay to an input to ff_s from f_e to make it like:

ff_s = (min(max((i_2 + w_5 * f_e(t-delay) + w_6 * f3(t) - w_7 * f_v), threshold), 25))

For now, as I understand correctly each solution is calculated directly from the previous timepoint solution (part where f_e, f_p, f_s, f_v = y).

Is it possible to somehow return a history of all previous timepoints or in other way introduce delay?

I have tried to search in the related topics at stackoverflow and did not find an answer.