Hi Daniel,
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
Your help would be highly appreciated.
Best,
Biswa
def taum(V):
from math 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 math 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 math 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 math 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 minfy
def hinfy(V):
from math import exp
ah = 0.07 * exp(-(V+65)/20)
bh = 1.0 / (1.0 + exp(-(V+35)/10.0))
hinit = ah/(ah + bh)
return hinfy
def ninfy(V):
from math 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 ninfy__author__="Biswa Sengupta"
__date__ ="$Jun 1, 2011 10:47:37 AM$"
if __name__ == "__main__":
from fipy import *
from scipy import *
from math import *
from numpy 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()
On 1 Jun 2011, at 16:55, Daniel Wheeler wrote:
>
> Biswa, You could try using fipy in conjunction with scipy's ode
> module. The benefit of this approach is that it will be failrly easy
> to set up and see how it works. You'll get a feeling for how useful
> the approach is in fairly short order. Cheers
>
> On Tue, May 31, 2011 at 12:14 PM, Biswa Sengupta <[email protected]> wrote:
>> Dear all,
>> I am a new user of FiPy and I thought it would be useful to ask the group
>> about my specific problem before I get into the nitty-gritty of learning
>> FiPy. Is it is possible to solve a coupled time-varying PDE along with a set
>> of state dependent ODEs, just like the Hodgkin-Huxley cable equations in
>> neuroscience (http://www.neuron.yale.edu/course/hhcable.html).
>> Has anyone in the group solved similarly coupled equations using FiPy? Also,
>> it would be great if there are any automatic differentiation routines that
>> could be applied to the data-structures used in FiPy?
>> Best,
>> Biswa
>
>
>
> --
> Daniel Wheeler
>
>
