Solving an ODE with bounds on a variable

Question

I am modelling a bacterial population that grows and dies while consuming a substrate. The model is defined by three ordinary differential equations, one for population, two for each substrate.

My issue is that I need to cap the maximum population to a certain integer. I have down this by adding a maximum_density_penalty to my growth rate, which maps the difference between the current population count and the maximum count to a float between 0 and 1 using the logistic function.

I also need all those values to always be greater than 0. I believe the equations to solve for this, but my issue, I believe, is that the t time step value given by the solve_ivp function is too large and overshoots my calculations. I can show full code for the model if necessary.

Below is a minimal reproducible example.

import numpy as np
from math import floor
from scipy.integrate import solve_ivp

# Initial conditions
initial_carbon_content = 1000000000000
initial_inorganic_carbon_content = 0
initial_cell_count = 10 ** 5

# Community paramaters
natural_death_fraction = 0.001 # %/day

# Heterotroph organism paramaters
ks = 882000000
carbon_content_per_cell = 100
growth_rate = 1  # day^-1 - Max growth rate.


# System paramaters
maximum_cell_count = 10 ** 9
duration = 11000 * 365.25

# Set up initial condition vector
y0 = [initial_carbon_content, initial_cell_count]
t_span = (0, duration)  # time grid

def model(t, y):
    organic_carbon_content = y[0]
    cell_count = y[1]

    # Growth
    x = max(maximum_cell_count - cell_count, 0)
    maximum_density_penalty = np.arctan(x)/(np.pi/2)
    penalized_growth_rate = growth_rate * maximum_density_penalty
    growth = floor(penalized_growth_rate * (organic_carbon_content / (ks + organic_carbon_content)) * cell_count)

    # Calculate the amount of cells at next step:
    # next_cell_count = (cell_count + growth - natural_deaths) * t

    # Natural Death
    natural_deaths = floor(natural_death_fraction * cell_count)

    # Organic carbon requirement
    required_organic_carbon_per_cell = 5
    required_organic_carbon = required_organic_carbon_per_cell * cell_count

    # Starvation
    starvation_deaths = floor(max((required_organic_carbon - organic_carbon_content), 0) / required_organic_carbon_per_cell)

    # Deaths
    deaths = starvation_deaths + natural_deaths

    # Net cell count change
    dndt = growth - deaths

    ## ORGANIC CARBON
    # Organic carbon outputs
    consumed_organic_carbon_from_population = required_organic_carbon_per_cell * (cell_count - starvation_deaths)
    consumed_organic_carbon_from_growth = dndt * carbon_content_per_cell
    total_organic_carbon_outputs = consumed_organic_carbon_from_population + consumed_organic_carbon_from_growth

    # Net organic carbon change
    dsdt = - total_organic_carbon_outputs

    assert 0 <= organic_carbon_content + dsdt
    assert 0 <= cell_count + dndt <= maximum_cell_count

    return [dsdt, dndt]

# solve the ODEs
soln = solve_ivp(fun=model, y0=y0, t_span=t_span)