Hello all,

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[2] * u[4] - u[2]
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[1]
maintenance_per_cell = p[2]
dOC_per_cell = p[3]
carrying_capacity = p[4]
pOC_input_rate = p[5]
dOC_input_rate = p[6]
inorganic_carbon_input_rate = p[7]
inorganic_carbon_fixing_rate = p[8]
inorganic_carbon_per_cell = p[9]
eea_rate = p[10]
Ks = p[11]
# Load state conditions
pOC_content = u[1]
dOC_content = u[2]
inorganic_carbon_content = u[3]
cell_count = u[4]
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[1] = pOC_input_rate - eea_rate*cell_count
du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_rate*cell_count
du[4] = 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