I am working on a mathematical model that describes the growth of 4 different bacterial populations and immune system cells under certain conditions. The model is governed by the equations below.

**POPULATION 1:**

$ {\underbrace{\frac{dN_{PS}}{dt}}_{\text{ Population change rate }}=\underbrace{r N_{PS}}_{\text{ Exponential growth }}\cdot\underbrace{\left(1-\frac{N_{PS}+N_{PR}}{K}\right)}_{\text{ Growth limitation }}-\underbrace{\theta_{PS}N_{PS}}_{\text{ Natural death }}-\underbrace{A_{PS}N_{PS}}_{\text{Biofilm formation}}+\underbrace{D_{BS}N_{BS}}_{\text{Biofilm dispersion}}-\underbrace{\phi_{PS}N_{PS}}_{\text{Mutation rate}}-\underbrace{\eta\delta_{PS}A_{m}N_{PS}}_{\text{Antibiotic inhibition}}-\left\{ \underbrace{\Gamma N_{PS}I}_{\text{Immune system}}\right\} }$

**POPULATION 2:**

$ {\frac{dN_{BS}}{dt}=(r-c_{b})N_{BS}\cdot\left(1-\frac{N_{BS}+N_{BR}}{K}\right)-\theta_{BS}N_{BS}+A_{PS}N_{PS}-D_{BS}N_{BS}-\phi_{BS}N_{BS}-\eta\delta_{BS}A_{m}N_{BS}-\left\{ \Gamma N_{BS}I\right\} }$

**POPULATION 3:**

$ {\frac{dN_{PR}}{dt}=(r-c_{r})N_{PR}\cdot\left(1-\frac{N_{PS}+N_{PR}}{K}\right)-\theta_{PR}N_{PR}-A_{PR}N_{PR}+D_{BR}N_{BR}+\phi_{PS}N_{PS}-\eta\delta_{PR}A_{m}N_{PR}-\left\{ \Gamma N_{PR}I\right\} }$

**POPULATION 4:**

$ {\frac{dN_{BR}}{dt}=(r-(c_{b}+c_{r}))N_{BR}\cdot\left(1-\frac{N_{BS}+N_{BR}}{K}\right)-\theta_{BR}N_{BR}+A_{PR}N_{PR}-D_{BR}N_{BR}+\phi_{BS}N_{BS}-\eta\delta_{BR}A_{m}N_{BR}-\left\{ \Gamma N_{BR}I\right\} }$

**NAIVE CELLS (IMMUNE SYSTEM):**

$ {\frac{dV}{dt}=\frac{-\sigma VB}{\pi+B}}$

**EFFECTOR CELLS (IMMUNE SYSTEM):**

$ {{\frac{dE}{dt}=(2\sigma V+\sigma E)\cdot\frac{B}{\pi+B}-hE\left(1-\frac{B}{\pi+B}\right)}}$

**MEMORY CELLS (IMMUNE SYSTEM):**

$ {{\frac{dM}{dt}=fEh\left(1-\frac{B}{\pi+B}\right)}}$

**TOTAL BACTERIAL DENSITY:**

$ {\frac{dB}{dt}=N_{PS}+N_{PR}+N_{BS}+N_{BR}}$

**IMMUNE SYSTEM DENSITY:**

$ {\frac{dI}{dt}=V+E+M}$

**ANTIBIOTIC UPTAKE 1:**

$ {{\eta}=\begin{cases} 1 & t_{1}\leq t\leq t_{1}+t_{2}\ 0 & t_{1}>t\:or\:t>t_{1}+t_{2} \end{cases}}$

**ANTIBIOTIC UPTAKE 2:**

$ {{\eta}=\begin{cases} 1 & B\geq\varOmega\ 0 & B<\varOmega \end{cases}}$

I am interested in how $ N_{PS}, N_{PR}, N_{BS}, N_{BR}, V, E, M $ change over time for which I implemented an algorithm in Python to solve the equations. Greek letters and other undescribed parameters (e.g. r, K, etc.) are mainly constants that are set at the beginning of the execution of the program.

One example of the functions used is shown below. Currently as you can also see from the code I’m using an Euler method to solve the equations, however I’d like to implement at least Heun’s method or even some higher order Runge-Kutta method to solve them.

I’m stuck already with Heun’s method and don’t know how to implement it. I ask for help with how to modify the following code and change it to for example Heun’s if that’s even possible with those equations.

`def sensitive_bacteria_PS(previous_density, growth_rate, density_PR, maximum_density_P, death_rate_PS, attachment_rate_PS, mutation_rate_PS, antibiotic_uptake, antibiotic_inhibition_PS, antibiotic_concentration, lymphocyte_inhibition, total_immune_cells, detachment_rate_BS, density_BS, time_step): return previous_density + ((previous_density * (growth_rate * (1 - (previous_density + density_PR) / maximum_density_P) - death_rate_PS - attachment_rate_PS - mutation_rate_PS - antibiotic_uptake * antibiotic_inhibition_PS * antibiotic_concentration - lymphocyte_inhibition * total_immune_cells) + detachment_rate_BS * density_BS) * time_step) `