I am trying to simulate cell growth, substrate consumption, and product in a fed-batch (a system where the substrate is added in different ways; exponential, continuous or bolus), I am using a system of equations based on monod ,I am using scipy to solve a system to simulate cell growth, i want to add substrate at certain time points, I tried to interpolate the value and add them to the function of the substrate, but the substrate is not being consumed

I tried interpolating the values of the substrate and adding them to the model, hoping the substrate would be consumed, but this is not happening. The substrate is accumulating, the cell concentrations do not increase, and the substrate is not consumed. I look some examples on GitHub, the difference is the feeding strategy they used doesn’t work for me Also, some examples on stack for loop iteration may solve the issue

I am expecting the simulation to be like this

import numpy as np
import pandas as pd
import scipy.integrate
from scipy.interpolate import interp1d
import plotly.graph_objects as go

umax = 0.0306
udmax = 0.02
Klys = 0.007
ks= 0.667
kipro= 3.6
Yx_s = 534000000
Yx_p = 5640000
Ypro_s = 0.0970
Cs = 500
F = 0
Fout = 0
V = 15
s= 5
Xv = 1502000000
Xd = 19000000
p= 0
pro = 0

x0 = [V, s, Xv, Xd, p, pro]  # State vector

# Time points and corresponding s feed rates
T_values = [0, 28.64, 100.57, 126.6, 158.6, 172.01, 200.64, 269.56, 298.57, 325]
Sf = [0, 4.856, 8.525, 6.984, 4.189, 6.367, 4.909, 7.761, 8.858, 0]

# Interpolate sF values over the time span
s_interpolator = interp1d(T_values, Sf, fill_value="extrapolate")

# bior function
def bior(t, y):
    V, s, Xv, Xd, pro, lac = y

    u = umax * (s / (ks + s)) * (kipro / (kipro + pro))
    ud = udmax * (pro / (pro + kipro))
    Qsubs = u / Yx_p
    Qp = u / Yx_p
    Qpro = Ylpro_p * Qp

    SF = s_interpolator(t)

    dV = F - Fout
    Dglu = -(Qs* Xv) + SF
    dXv = (u - ud) * Xv
    DXd = ud * Xv - Klys * Xd
    Dp = Qp * Xv
    Dpro= Qpro * Xv

    return [dV, Dgs, dXv, DXd, Dmp, Dpro]

# Time span
tspan = (0, 325)
tsmooth = np.linspace(0, 325, 1000)

sol = scipy.integrate.solve_ivp(bior, tspan, x0, t_eval=tsmooth)
# Create DataFrame
df = pd.DataFrame({
    't': sol.t,
    'V': sol.y[0],
    's': sol.y[1],
    'Xv': sol.y[2],
    'Xd': sol.y[3],
    'mab': sol.y[4],
    'lac': sol.y[5]
})


fig1 = go.Figure()


fig1.add_trace(go.Scatter(x=df['t'], y=df['Xv'], mode='lines', name='Xv'))

fig1.update_layout(
    title='X Concentration over Time',
    xaxis_title='Time',
    yaxis_title='X Concentration'
)


fig1.show()

fig2 = go.Figure()


fig2.add_trace(go.Scatter(x=df['t'], y=df['s'], mode='lines', name='s'))

# Update layout for glucose plot
fig2.update_layout(
    title='s Concentration over Time',
    xaxis_title='Time',
    yaxis_title='Glucose Concentration'
)

fig2.show()

print(df)