优化功能参数

Question

我简单解释一下附加的程序代码应该做什么。我们在 runs = 100 之前给出了一些传球。我们给出 I = 10.

例如我们设置area_factor = 1。然后函数 HH_model(I,area_factor) 执行以下操作：运行 100 次这个 I 和这个 area_factor 和 return 屏障 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))

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)

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

Answer 1

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

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

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

优化功能参数

Optimizing function parameters

python

modeling