I explain briefly what the attached program code should do. We give a number of passes before runs = 100. and we give I = 10.
For example we set the area_factor = 1. Then the function HH_model(I,area_factor) does the following:
run 100 times with this I and this area_factor and return the number of times the barrier 60 is broken -- this is checked in the if max(v[:]-v_Rest) > 60 query.
Now I want to do the following: Determine that area_factor so that the number of count matches observations as well as possible. For example, I know from measurements
HH_model(2*I,area_factor) = 70
HH_model(I,area_factor)=50
HH_model(0.5*I,area_factor) = 30
...
how can I find the area_factor for a given I, so that the difference to the observations becomes minimal.
import matplotlib.pyplot as py
import numpy as np
import scipy.optimize as optimize
# HH parameters
v_Rest = -65 # in mV
gNa = 120 # in mS/cm^2
gK = 36 # in mS/cm^2
gL = 0.3 # in mS/cm^2
vNa = 115 # in mV
vK = -12 # in mV
vL = 10.6 # in mV
#Number of runs
runs = 30
c = 1 # in uF/cm^2
#performing bisection-procedure
ROOT = True
def HH_model(I,area_factor):
count = 0
t_end = 10 # in ms
delay = 1 # in ms
duration = 0.3 # in ms
dt = 0.01 # in ms
I = I
area_factor = area_factor
#geometry
d = 2 # diameter in um
r = d/2 # Radius in um
l = 10 # Length of the compartment in um
A = (2 * np.pi * r * l * 1e-8)*area_factor # surface [cm^2]
C = c * A # uF
for j in range(0,runs):
# Introduction of equations and channels
def alphaM(v): return 12 * ((2.5 - 0.1 * (v)) / (np.exp(2.5 - 0.1 * (v)) - 1))
def betaM(v): return 12 * (4 * np.exp(-(v) / 18))
def betaH(v): return 12 * (1 / (np.exp(3 - 0.1 * (v)) + 1))
def alphaH(v): return 12 * (0.07 * np.exp(-(v) / 20))
def alphaN(v): return 12 * ((1 - 0.1 * (v)) / (10 * (np.exp(1 - 0.1 * (v)) - 1)))
def betaN(v): return 12 * (0.125 * np.exp(-(v) / 80))
# compute the timesteps
t_steps= t_end/dt+1
# Compute the initial values
v0 = 0
m0 = alphaM(v0)/(alphaM(v0)+betaM(v0))
h0 = alphaH(v0)/(alphaH(v0)+betaH(v0))
n0 = alphaN(v0)/(alphaN(v0)+betaN(v0))
# Allocate memory for v, m, h, n
v = np.zeros((int(t_steps), 1))
m = np.zeros((int(t_steps), 1))
h = np.zeros((int(t_steps), 1))
n = np.zeros((int(t_steps), 1))
# Set Initial values
v[:, 0] = v0
m[:, 0] = m0
h[:, 0] = h0
n[:, 0] = n0
### Noise component
knoise= 0.003 #uA/(mS)^1/2
### --------- Step3: SOLVE
for i in range(0, int(t_steps)-1, 1):
# Get current states
vT = v[i]
mT = m[i]
hT = h[i]
nT = n[i]
# Stimulus current
IStim = 0
if delay / dt <= i <= (delay + duration) / dt:
IStim = I * A # in uA
else:
IStim = 0
# Compute change of m, h and n
m[i + 1] = (mT + dt * alphaM(vT)) / (1 + dt * (alphaM(vT) + betaM(vT)))
h[i + 1] = (hT + dt * alphaH(vT)) / (1 + dt * (alphaH(vT) + betaH(vT)))
n[i + 1] = (nT + dt * alphaN(vT)) / (1 + dt * (alphaN(vT) + betaN(vT)))
# Ionic currents
iNa = gNa * m[i + 1] ** 3. * h[i + 1] * (vT - vNa)
iK = gK * n[i + 1] ** 4. * (vT - vK)
iL = gL * (vT-vL)
Inoise = (np.random.normal(0, 1) * knoise * np.sqrt(gNa * A))
IIon = ((iNa + iK + iL) * A) + Inoise #
# Compute change of voltage
v[i + 1] = vT + ((-IIon + IStim) / C) * dt # in ((uA / cm ^ 2) / (uF / cm ^ 2)) * ms == mV
# adjust the voltage to the resting potential
v = v + v_Rest
# test if there was a spike
if max(v[:]-v_Rest) > 60:
count += 1
return count
Ich habe folgendes versucht:
I = 30
xdata = np.array([0.92*I,I,1.05*I])
ydata = np.array([28,100,110])
y0=np.array([1,1,1])
def g(y,xdata,ydata):
return ydata - HH_model(xdata,y)
fit = optimize.leastsq(g, y0, args=(xdata, ydata))
File "", line 126, in HH_model v[i + 1] = vT + ((-IIon + IStim) / C) * dt
ValueError: could not broadcast input array from shape (3) into shape (1)
how can I get around this and make the input in the correct format?
The result of your line 126 is a three dimensional array with three times the same value. This size-3 array does not fit into an element of v, which has size-1 elements as you initialized them this way.
Therefore, you could add a [0]:
v[i + 1] = (vT + ((-IIon + IStim) / C) * dt)[0]
Furthermore, I think you do not need to allocate memory. You could for example use numpy.append in line 126.