优化函数参数

Question

我简要解释了附加的程序代码应该做什么。我们在runs = 100 之前给出了一些通行证。我们给I = 10。

例如我们设置area_factor = 1。然后函数HH_model(I,area_factor) 执行以下操作： 使用这个 I 和这个 area_factor 运行 100 次，并返回障碍 60 被打破的次数——这在 if max(v[:]-v_Rest) > 60 查询中检查。

现在我想做以下事情：确定 area_factor 以便计数的数量尽可能匹配观察结果。 例如，我从测量中知道

HH_model(2*I,area_factor) = 70
HH_model(I,area_factor)=50
HH_model(0.5*I,area_factor) = 30

...

如何找到给定 I 的 area_factor，从而使观察结果的差异变得最小。

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))

文件“”，行 126，在 HH_model v[i + 1] = vT + ((-IIon + IStim) / C) * dt

ValueError: 无法将输入数组从形状 (3) 广播到形状 (1)

我怎样才能解决这个问题并以正确的格式输入？

Answer 1

第 126 行的结果是一个具有三倍相同值的三维数组。这个大小为 3 的数组不适合 v 的元素，当您以这种方式初始化它们时，它具有大小为 1 的元素。

因此，您可以添加 [0]： v[i + 1] = (vT + ((-IIon + IStim) / C) * dt)[0]

此外，我认为您不需要分配内存。例如，您可以在第 126 行使用 numpy.append。