On Jun 2, 2011, at 5:41 AM, Biswa Sengupta wrote:
> As per your advice I tried to pose my set of PDE-ODE system as a PDE + 3
> degenerate PDEs (i.e., the ODEs without the space term). I have tried to work
> along the examples for time marching but certainly I am not getting something
> right. I attach herewith the code that I am trying to use. Can you spot some
> obvious mistake. I get the following error:
>
> 4565:src biswasengupta$ python cable.py
> Traceback (most recent call last):
> File "cable.py", line 63, in <module>
> mEq = TransientTerm(coeff=taum(V)) == ImplicitSourceTerm(minfy(V) - m)
> File "/Users/biswasengupta/Desktop/finelPDE/python/cable/src/actinact.py",
> line 3, in taum
> am = 0.1*(V + 40) / (1-exp(-(V+40)/10.0))
> File
> "/Library/Frameworks/EPD64.framework/Versions/7.0/lib/python2.7/site-packages/FiPy-2.1.2-py2.7.egg/fipy/variables/variable.py",
> line 1146, in __float__
> return float(self.getValue())
> TypeError: only length-1 arrays can be converted to Python scalars
In your actinact.py, replace all instances of
from math import exp
with
from fipy.tools.numerix import exp
The Python math module only understands single floating point numbers and
should never be used with FiPy. FiPy's numerix, based on NumPy, can deal with
arrays of floats.
----
In cable.py, you should remove
from math import *
for the same reason. Also
from numpy import *
is not necessary because all of the same code (with appropriate modifications)
is obtained when you import fipy. By importing numpy directly, you override
some changes that FiPy requires. Finally,
from scipy import *
doesn't do anything. You must import the specific submodules you want from
SciPy.
----
Once you're past that, you will encounter a few more issues before your script
will run, some of them quite obscure:
* TypeError: unsupported operand type(s) for *: 'function' and 'float'
This turns out to be caused by
return ninfy
on actinact.py line 41. Change this to
return ninit
and likewise for minfy and hinfy.
* self.matrix.matvec(other, y)
ValueError: arg 1 must be a 1-dimensional double array of appropriate size.
This is because V holds an array of integers and it needs to be floats.
Change cable.py
line 42 to
restvol = -65.
Probably a good idea to change valueLeft and valueRight (but not absolutely
necessary).
* Finally, the problem runs, but doesn't show anything. This is because
axonlen = 1
dx = axonlen/nx
gives a grid spacing of zero because of integer math. Change to
axonlen = 1.
The problem now runs. Whether it gives sensible answers is up to you to figure
out! 8^)
I've attached my modified scripts.
__author__="Biswa Sengupta"
__date__ ="$Jun 1, 2011 10:47:37 AM$"
if __name__ == "__main__":
from fipy import *
from actinact import *
###########################################################################
# squid axon constants
###########################################################################
PI = 3.141528 # dimensionless
cable_length = 1 # in cm
cable_radius = 0.159155 # in cm
patch_area = (2*PI*cable_radius*cable_length) # in cm^2
C_m_spec = 1 # uF/cm^2
cm = C_m_spec*patch_area # uF
e_Na = 50 # in mV
gNa_spec = 120 # in mS/cm^2
gNa = gNa_spec*patch_area # in mS
e_K = -77 # in mV
gK_spec = 20 # in mS/cm^2
gK = gK_spec*patch_area # in mS
e_l = -54.4 # in mV
g_l_spec = 0.3 # in mS/cm^2
gL = g_l_spec*patch_area # in mS
idc_spec = 20 # in uA/cm^2
idc = idc_spec*patch_area # in uA/cm^2
ra_spec = 35.4 # in Ohm cm
ra = cable_length*(ra_spec/(PI*cable_radius**2)) # in Ohm
lamda = sqrt((cable_radius)/(2.0*ra_spec*g_l_spec*(10**-3))) # cm
taumem = cm/gL # ms
###########################################################################
nx = 100
axonlen = 1.
dx = axonlen/nx
mesh = Grid1D(nx = nx, dx = dx)
restvol = -65.
V = CellVariable(name="Voltage", mesh=mesh, value=restvol,hasOld=1)
m = CellVariable(name="m", mesh=mesh, value=minfy(restvol),hasOld=1)
h = CellVariable(name="h", mesh=mesh, value=hinfy(restvol),hasOld=1)
n = CellVariable(name="n", mesh=mesh, value=ninfy(restvol),hasOld=1)
valueLeft = -65
valueRight = -65
BCs = (FixedValue(faces=mesh.getFacesRight(), value=valueRight),FixedValue(faces=mesh.getFacesLeft(), value=valueLeft))
# ionic currents
ina = gNa*m*m*m*h*(V-e_Na)
ik = gK*n*n*n*n*(V-e_K)
il = gL*(V-e_l)
ionic = ina + ik + il
# voltage evolution
volEq = TransientTerm(coeff=taumem) == ImplicitDiffusionTerm(coeff=lamda**2) - ImplicitSourceTerm(coeff=ionic)
# the (in) activation ODEs described as degenerate PDEs
mEq = TransientTerm(coeff=taum(V)) == ImplicitSourceTerm(minfy(V) - m)
hEq = TransientTerm(coeff=tauh(V)) == ImplicitSourceTerm(hinfy(V) - h)
nEq = TransientTerm(coeff=taun(V)) == ImplicitSourceTerm(ninfy(V) - n)
timeStepDuration = .001
steps = 100
viewer = Viewer(vars=V)
for i in range(steps):
V.updateOld()
m.updateOld()
h.updateOld()
n.updateOld()
volEq.sweep(V,dt=timeStepDuration)
mEq.sweep(m,dt=timeStepDuration)
hEq.sweep(h,dt=timeStepDuration)
nEq.sweep(n,dt=timeStepDuration)
viewer.plot()def taum(V):
from fipy.tools.numerix import exp
am = 0.1*(V + 40) / (1-exp(-(V+40)/10.0))
bm = 4.0 * exp(-(V+65)/18.0)
taum = 1.0/(am + bm)
return taum
def tauh(V):
from fipy.tools.numerix import exp
ah = 0.07 * exp(-(V+65)/20)
bh = 1.0 / (1.0 + exp(-(V+35)/10.0))
tauh = 1.0/(ah + bh)
return tauh
def taun(V):
from fipy.tools.numerix import exp
an = 0.01*(V+55) / (1.0 - exp(-(V+55)/10))
bn = 0.125 * exp((-V-65)/80.0)
taun = 1.0/(an + bn)
return taun
def minfy(V):
from fipy.tools.numerix import exp
am = 0.1*(V + 40) / (1-exp(-(V+40)/10.0))
bm = 4.0 * exp(-(V+65)/18.0)
minit = am/(am + bm)
return minit
def hinfy(V):
from fipy.tools.numerix import exp
ah = 0.07 * exp(-(V+65)/20)
bh = 1.0 / (1.0 + exp(-(V+35)/10.0))
hinit = ah/(ah + bh)
return hinit
def ninfy(V):
from fipy.tools.numerix import exp
an = 0.01*(V+55) / (1.0 - exp(-(V+55)/10))
bn = 0.125 * exp((-V-65)/80.0)
ninit = an/(an + bn)
return ninit
