I have a set of relatively simple differential equations modeling growth of bacteria as a function of available substrates. Unfortunately, my model fails to be solved in one particular case, and I cannot figure out why. Here’s is the (simplified) model:
function solve_model(p, u0) p = [2.767317774014066e-7, 0.7433220714328409, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0] # FAILS # p = [0.06, 0.1806273414373398, 15.7, 1.0e9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.012192, 882000.0] # WORKS u0 = [1.6382099999999998e13, 3.4061512877350002e10, 0.0, 100000.0, 1.461e7] # Callback for the max() function in the model - causes a discontinuity. # If max equation changes, this condition will have to change. condition(u, t, integrator) = p * u - u affect!(integrator) = nothing max_cb = ContinuousCallback(condition, affect!, save_positions=(false,false)) # Build the ODE Problem and Solve prob = ODEProblem(model, u0[1:end-1], (0.0, last(u0)), p) sol = solve(prob, Rosenbrock23(), callback=max_cb, tstops=stops, isoutofdomain=is_invalid_domain, maxiters=1e7) return sol end # Implements the differential equations that define the model function model(du,u,p,t) # Paramaters mu_max = p maintenance_per_cell = p dOC_per_cell = p carrying_capacity = p pOC_input_rate = p dOC_input_rate = p inorganic_carbon_input_rate = p inorganic_carbon_fixing_rate = p inorganic_carbon_per_cell = p eea_rate = p Ks = p # Load state conditions pOC_content = u dOC_content = u inorganic_carbon_content = u cell_count = u growth = mu_max * (dOC_content / (Ks + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity) required_dOC = maintenance_per_cell * cell_count dOC_missing = max(required_dOC - dOC_content, 0) starvation_deaths = dOC_missing == 0 ? 0 : dOC_missing / required_dOC_per_cell deaths = starvation_deaths dOC_consumption = maintenance_per_cell * (cell_count - deaths) eea_rate = eea_rate * (1 / (1 + exp(-pOC_content))) * 2 - eea_rate du = pOC_input_rate - eea_rate*cell_count du = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count du = growth - deaths end
I originally believed that the reason instability caused the solver to abort was due to discontinuities at the max function. However implementing fixes for that did not solve the issue, see here for attempted solutions. This includes the recommended first steps of trying a different solver and lowering tolerance, the issue is definitely my model.
I was hoping more experienced modelers than I on this forum could offer some direction as to why this set of parameter fails. And from there, I hope to resolve many runs aborting when I conduct a sensitivity analysis. One thing I am aware of is that the model never seems to settle on zero values. Often it gives approximations such as 1e-300.
I would also appreciate general debugging tips, as stepping through 1000 iterations of a model call is impractical and so is manual logging. What are some strategies used?
Full model here