dynamic flux balance analysis of genome-scale metabolic models in COBRApy

Question

I am trying to simulate the metaboic behavior of GEMs over a specific period of time using dynanuc flux balance analysis. But I am running into some problems. First of all I am not using the dfba package directly for the simulation, I am using ordinary differential equations to solve the FBA values of the model at each time step, but after several attempts I found that the run is very slow, and after trying to shorten the time intervals and reduce the time step sizes I found that the substrate glucose concentration went nagetiva, which is obviously not feasible, and I suspect that it may be due to my incorrect setting of the parameters, but maybe it's some other problem, hopefully someone can help me with it.

Here is my code:

#read model
model_sa = cobra.io.read_sbml_model('/Users/lishijia/Downloads/iYS854.xml')

#set initial condition
co_parameters = {'sa_h_k':5,'sa_h_y':1,
                 'sa_nh4_k':5,'sa_nh4_y':1,
                 'sa_o2_k':5,'sa_o2_y':1,
                 'sa_glc_k':1,'sa_glc_y':1,
                 'sa_co2_y':-1,'sa_glyclt_y':-1,
                 'flow_rate':0.05}
#set the initial concentration of substrates, products, and biomass
initial_states = np.array([10,10,60,60,0,0,0.1])

#define growth rate function
def muSa(concentration_vec):

    consumable_glc = (concentration_vec[3]/concentration_vec[3]+co_parameters['sa_glc_k'])*10
    consumable_h = (concentration_vec[0]/(concentration_vec[0]+co_parameters['sa_h_k']))*5
    consumable_nh4 = (concentration_vec[1]/(concentration_vec[1]+co_parameters['sa_nh4_k']))*5
    consumable_o2 = (concentration_vec[2]/(concentration_vec[2]+co_parameters['sa_o2_k']))*10

    model_sa.reactions.EX_glc__D_e.lower_bound = -1*consumable_glc
    model_sa.reactions.EX_h_e.lower_bound = -1*consumable_h
    model_sa.reactions.EX_nh4_e.lower_bound = -1*consumable_nh4
    model_sa.reactions.EX_o2_e.lower_bound = -1*consumable_o2

    sol_sa = model_sa.optimize()

    consumed_glc = max(0,model_sa.reactions.EX_glc__D_e.flux*-1)
    consumed_h = max(0,model_sa.reactions.EX_h_e.flux*-1)
    consumed_nh4 = max(0,model_sa.reactions.EX_nh4_e.flux*-1)
    consumed_o2 = max(0,model_sa.reactions.EX_o2_e.flux*-1)

    produced_co2 = max(0,model_sa.reactions.EX_co2_e.flux)
    produced_glyclt = max(0,model_sa.reactions.EX_glyclt_e.flux)

    return (consumed_h,consumed_nh4,consumed_o2,consumed_glc,produced_co2,produced_glyclt,
            sol_sa.objective_value)

#solve model by ode
def ode(t,states):
    mu_sa = muSa(states)

    pop_sa = states[6]

    dSa = pop_sa*(mu_sa[-1]-co_parameters['flow_rate'])

    sa_h = -(mu_sa[0]/co_parameters['sa_h_y'])*pop_sa
    sa_nh4 = -(mu_sa[1]/co_parameters['sa_nh4_y'])*pop_sa
    sa_o2 = -(mu_sa[2]/co_parameters['sa_o2_y'])*pop_sa
    sa_glc = -(mu_sa[3]/co_parameters['sa_glc_y'])*pop_sa

    sa_co2 = -(mu_sa[4]/co_parameters['sa_co2_y'])*pop_sa
    sa_glyclt = -(mu_sa[5]/co_parameters['sa_glyclt_y'])*pop_sa

    derivatives = np.zeros(len(states))

    derivatives[0] = co_parameters['flow_rate']*(initial_states[0]-states[0])+sa_h
    derivatives[1] = co_parameters['flow_rate']*(initial_states[1]-states[1])+sa_nh4
    derivatives[2] = co_parameters['flow_rate']*(initial_states[2]-states[2])+sa_o2
    derivatives[3] = co_parameters['flow_rate']*(initial_states[3]-states[3])+sa_glc
    derivatives[4] = co_parameters['flow_rate']*(initial_states[4]-states[4])+sa_co2
    derivatives[5] = co_parameters['flow_rate']*(initial_states[5]-states[5])+sa_glyclt

    derivatives[6] = dSa

    if ode.pbar is not None:
        ode.pbar.update(1)
        ode.pbar.set_description('t = {:.9f}'.format(t))

    return derivatives

with tqdm() as pbar:
    ode.pbar = pbar
    solution3 = solver(fun=ode, t_span = (0, 10), y0 = initial_states, t_eval = np.linspace(0, 10, 100))