I've been thinking about moving my modelling from Mathematica (expensive, poorly managed activation, expensive) to Python/Scipy (free, friendly community, many scientific modules). In theory, all of the tools that I employ in Mathematica are available through the Scipy and Numpy Python modules.
I enjoy using Python for a number of reasons:
- It's free
- I am able to use TextMate, a fantastic text editor available on the Mac avoiding the clunky notebook entry of Mathematica
- The community is active and friendly
- Many, many tools are available for scientific research e.g. the excellent SciPy, and Biopython to name two that are relevant to my work
Granted that the visualisation tools aren't quite as user-friendly as those available in Mathematica and the symbolic computing power of Mathematica is sometimes missed, the tools available are more than adequate for my needs.
To directly compare the two systems, I obtained a simple Hodgkin-Huxley model from the Mathematica demonstration site titled: Neural Impulse: The Action Potential in Action. I stripped it of its interactive elements and plotting functions to the bare bones:
Vmi = 0; Gna = 120; Gk = 36; Gl = 0.3; Ena = 115; Ek = -12; El = 10.6;
alpham[Vm_] := (0.1 (25 - Vm))/(Exp[(25 - Vm)/10] - 1)
betam[Vm_] := 4 Exp[-Vm/18]
alphan[Vm_] := (0.01 (10 - Vm))/(Exp[(10 - Vm)/10] - 1)
betan[Vm_] := 0.125 Exp[-Vm/80]
alphah[Vm_] := 0.07 Exp[-Vm/20]
betah[Vm_] := 1/(Exp[(30 - Vm)/10] + 1)
gna[t_] := (Gna) ((m[t])^3) h[t]
gk[t_] := (Gk) (n[t])^4
duration = 20;
strength = 3
Istim[t_, duration_, strength_] := If[0 < t < duration, strength, 0]
For[i = 1, i <= 500, i ++,
NDSolve[{Vm'[t] ==
Istim[t, duration, strength] - 1 (gna[t] (Vm[t] - Ena)) -
1 (gk[t] (Vm[t] - Ek)) - (Gl (Vm[t] - El)),
m'[t] == (alpham[Vm[t]] (1 - m[t])) - (betam[Vm[t]] m[t]),
n'[t] == (alphan[Vm[t]] (1 - n[t])) - (betan[Vm[t]] n[t]),
h'[t] == (alphah[Vm[t]] (1 - h[t])) - (betah[Vm[t]] h[t]), Vm[0] == Vmi,
m[0] == alpham[Vmi]/(alpham[Vmi] + betam[Vmi]),
n[0] == alphan[Vmi]/(alphan[Vmi] + betan[Vmi]),
h[0] == alphah[Vmi]/(alphah[Vmi] + betah[Vmi])}, {Vm, m, n, h}, {t, 0,
50}]]
Time to completion: 9 seconds
Then I coded the same model in Python (source available here):
import pylab
import numpy
from scipy.integrate import odeint
import math
from numpy import array,linspace
G_L = 0.3
E_L = 10.6
G_Na = 120
E_Na = 115
G_K = 36
E_K = -12
e = math.e # Euler's constant
v0 = 0.0
alpham = (0.1*(25-v0)) / (e**((25-v0)/10)-1)
betam = 4 * e**(-v0/18)
alphah = 0.07 * e**(-v0 / 20)
betah = 1 / (e**((30-v0)/10)+1)
alphan = (0.01 * (10 - v0)) / (e**((10 - v0) / 10)-1)
betan = 0.125 * e**(-v0/80)
m0 = alpham/(alpham + betam)
n0 = alphan/(alphan + betan)
h0 = alphah/(alphah + betah)
def i_stim(t):
if t < 20: return 4
else: return 0.0
def alpha_m(x):
return (0.1 *(25-x))/(e**((25-x)/10)-1)
def beta_m(x):
return 4*e**(-x/18)
def alpha_h(x):
return 0.07 * e**(-x / 20)
def beta_h(x):
return 1.0 / (e**((30-x)/10)+1)
def alpha_n(x):
return (0.01 * (10-x)) / (e**((10-x)/10)-1)
def beta_n(x):
return 0.125 * e**(-x/80)
def diff_eqns(x,t):
"""Differential equations for the HH model"""
dm_dt = alpha_m(x[3]) * (1.0 - x[0]) - beta_m(x[3]) * x[0]
dh_dt = alpha_h(x[3]) * (1.0 - x[1]) - beta_h(x[3]) * x[1]
dn_dt = alpha_n(x[3]) * (1.0 - x[2]) - beta_n(x[3]) * x[2]
dv_dt = i_stim(t) - (G_L * (x[3] - E_L)) - (x[0]**3 * x[1] * G_Na * (x[3] - E_Na)) - (x[2]**4 * G_K * (x[3] - E_K))
return array([dm_dt, dh_dt, dn_dt, dv_dt])
t = linspace(0,50,1000)
i=1
while i < 500:
odeint(diff_eqns, array([m0,h0,n0,v0]), t)
i = i+1
Time to completion: 223 s
That is significantly slower than equivalent code in Mathematica. I enjoyed writing it in Python, but it was not difficult in Mathematica and considering the 10+ x speed up it does make me wonder whether I should make the switch. I have read a little about speeding up numerical solving in Python with the use of Fortran and C. But considering everything else that I have to learn it might just be better to stick with the tried and true. What do others think? Is there a simple way to speed up the Python code? I'll update the benchmarks as the code improves.
(Mac 10.6.2, Python 2.6.4 with 64 bit Numpy, SciPy; Mathematica 7.01)
Update: Now that I am running Mathematica 8, I thought I would repeat this very inaccurate benchmarking: time to completion 7s in Mathematica 8. A modest, but welcome improvement.
If you are new to Mathematica, the best help I have had after the included documentation is Mathematica Cookbook
by Sal Mangano; well worth the small price.
