Using Sympy for Biological ODEs Exponential and Logistic Growth solutions

I’ve been getting more into theoretical/analytical modeling lately, and I’ve been playing around with software to help me do some of the more complicated calculus involved (for two reasons: 1) I’m lazy, 2) My skills are very rusty). Python (of course) provides excellent symbolic capabilities through the Sympy module. I played around with it to try and get a feel for how it works, but I couldn’t find any help online, nor anyone who has posted a tutorial, on analyzing basic biological/ecological models with Sympy. So here is my version. Below, I solve both the exponential growth and logistic growth models using Sympy, then plot the results. Here is a step-by-step tutorial for Bio-Sympy.

Exponential Growth

Exponential growth is defined by the differential equation

\( \frac{dy(t)}{dt} = k*y(t) \)

and this ODE has the analytical solution of

\( y(t) = Y_0 e^{kt} \)

So how do we use Sympy to go from the ODE to the general solution? Well, like this.

First, gotta import the modules we need, and then initialize pretty printing just so the output is readable (see the Sympy docs for this):

import sympy as sm
import numpy as np
import matplotlib.pyplot as plt
sm.init_printing()

Next, define the variables k (the intrinsic growth rate) and t (for time). Then, make y a function of t

from sympy.abc import t, k
y = sm.Function('y')(t)

Then, define the derivative of y with respect to t (the left-hand side of the ODE), and then define the right-hand side of the ODE.

dy = y.diff(t)
rhs = k*y

We need to set these two quantities equal to one another, as in the ODE above, which we can do using a Sympy Equality

eq = sm.Eq(dy, rhs)

If you print eq, it should give you the differential equation. Now that we have the differential equation set up, all we need to do is solve. Since this is a simple one, Sympy can do it on its own with no hints or guesses:

sol = sm.dsolve(eq)
sol

Where sol should give you the analytical solution \( y(t) = C_1 e^{kt} \). But \( C_1 \) here is just a constant, we want to put it in terms of the initial conditions. (I know this is trivial in this case, but bear with me). First, we need to find the initial conditions. We do so by substituting 0 in for t, and then setting that equal to n0.

t0 = sol.args[1].subs({'t': 0})
n0 = sm.symbols('n0')
eq_init = sm.Eq(n0, t0)

The object eq_init should now be \( n_0 = C_1 \). Great. Now, solve that in terms of \( C_1 \) and substitute it back into the original equation.

C1 = eq_init.args[1] # this chunk just isolates the C1 variable as a symbol so we can use it in Sympy
sm.solve(eq_init, C1)
init_solve= sm.solve(eq_init, C1)
final = sol.subs(C1, init_solve[0])

final is now \(  y(t) = n_0 e^{kt} \). Just like it should be! I’ll leave plotting that one up to you.

Logistic Growth

Logistic growth is defined by the differential equation

\( \frac{dy(t)}{dt} = ry(t)\left ( 1 -\frac{y(t)}{K} \right ) \)

The steps for solving this are identical to the steps above:

# define the needed parameters
from sympy.abc import r, K, t
y = sm.Function('y')(t)
# define the differential equation
dy = y.diff(t)
rhs = r*y*(1 - y/K)
# define the equality
eq = sm.Eq(dy, rhs)
# find the general solution
sol = sm.dsolve(eq))

sol is now more complex:

\( y(t) = \frac{Ke^{C_1 K + rt}}{e^{C_1 K + rt} – 1} \)

# find out what C is in terms of y0
t0 = sol.args[1].subs({'t':0})
n0 = sm.symbols('n0')
eq_init = sm.Eq(n0, t0)
# it takes a little more work to isolate C1 here
# try each step for yourself to see what it does
C1 = t0.args[2].args[0].args[0]
t0_sol = sm.solve(eq_init, C1)

The initial conditions are also rather complex:

\( C_1 =\frac{\log \left (  \frac{-n_0)}{K-n_0} \right )}{K} \)

But we can substitute that back into the original solution to get the general solution in terms of the initial conditions

# substitute the expression for C1 back into the equation
final = sol.args[1].subs(C1, t0_sol[0])

where final is now

\( y(t) = \frac{-K n_0 e^{rt}}{  (K-n_0) \left ( – \frac{n_0 e^{rt}}{K-n_0} – 1  \right ) } \)

Holy cow! That’s a doozy, and it doesn’t look anything like the general solutions you’ll find elsewhere. Easy fix, just ask Sympy to simplify it!

final.simplify()

and we have

\( y(t) = \frac{ K n_0 e^{rt} }{  K + n_0 e^{rt} – n_0  } \)

which we can easily simplify further in the denominator ourselves:

\( y(t) = \frac{ K n_0 e^{rt} }{  K + n_0 (e^{rt} -1) } \)

and we have the general solution to logistic growth that everyone knows and love!

Viola, Sympy for solving biological ODEs