Grid Search over function

Question

The function HH_model(I,area_factor) has as return value the number of spikes which are triggered by n runs. Assuming 1000 runs, there are 157 times that max(v[]-v_rest) > 60, then the return value of HH_model(I,area_factor) is 157. Now I know value pairs from another model - the x-values are related to the stimulus I, while the y-values are the number of spikes.

I have written these values as a comment under the code. I want to choose my input parameters I and area_factor in a way that the error to the data is as small as possible. I have no idea how I should do this optimization.

import matplotlib.pyplot as py
import numpy as np
import scipy.optimize as optimize

# HH parameters
v_Rest = -65    # in mV
gNa = 1200      # in mS/cm^2
gK = 360      # in mS/cm^2
gL = 0.3*10      # in mS/cm^2
vNa = 115      # in mV
vK = -12       # in mV
vL = 10.6      # in mV

#Number of runs

runs = 1000


c = 1         # in uF/cm^2

ROOT = False


def HH_model(I,area_factor):
    
    count = 0
    t_end = 10  # in ms
    delay = 0.1     # in ms
    duration = 0.1    # in ms
    dt = 0.0025   # in ms
    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 = (1*10**(-8))*area_factor      # surface [cm^2]
    I = I
    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.0005  #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    # 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)[0]   # 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


# some datapoints from another model out of 1000 runs. ydata means therefore 'count' out of 1000 runs. 
# xdata = np.array([0.92*I,0.925*I,0.9535*I,0.975*I,0.9789*I,I,1.02*I,1.043*I,1.06*I,1.078*I,1.09*I])
# ydata = np.array([150,170,269,360,377,500,583,690,761,827,840])

EDIT:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize

# 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 = 1000


c = 1         # in uF/cm^2



def HH_model(x,I,area_factor):
    
    count = 0
    t_end = 10  # in ms
    delay = 0.1     # in ms
    duration = 0.1    # in ms
    dt = 0.0025   # in ms
    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 = (1*10**(-8))*area_factor      # surface [cm^2]
    I = I*x
    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.0005  #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    # 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)[0]   # 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


def loss(parameters, model, x_ref, y_ref):
    # unpack multiple parameters
    I, area_factor = parameters
    # compute prediction
    y_predicted = np.array([model(x, I, area_factor) for x in x_ref])
    # compute error and use it as loss
    mse = ((y_ref - y_predicted) ** 2).mean()
    return mse


# some datapoints from another model out of 1000 runs. ydata means therefore 'count' out of 1000 runs. 
xdata = np.array([0.92,0.925,0.9535, 0.975, 0.9789, 1])
ydata = np.array([150,170,269, 360, 377, 500])

y_data_scaled = ydata / runs

y_predicted = np.array([HH_model(x,I=10**(-3), area_factor=1) for x in xdata])

parameters = (10**(-3), 1)
mse0 = loss(parameters, HH_model, xdata, y_data_scaled)


# compute the parameters that minimize the loss (alias, the error between the data and the predictions of the model)
optimum = minimize(loss, x0=np.array([10**(-3), 1]), args=(HH_model, xdata, y_data_scaled))

# compute the predictions with the optimized parameters
I = optimum['x'][0]
area_factor = optimum['x'][1]
y_predicted_opt = np.array([HH_model(x, I, area_factor) for x in xdata])

# plot the raw data, the model with handcrafted guess and the model with optimized parameters
fig, ax = plt.subplots(1, 1)
ax.set_xlabel('input')
ax.set_ylabel('output predictions')
ax.plot(xdata, y_data_scaled, marker='o')
ax.plot(xdata, y_predicted, marker='*')
ax.plot(xdata, y_predicted_opt, marker='v')
ax.legend([
    "raw data points",
    "initial guess",
    "predictions with optimized parameters"
])

user_14637936 · Answer

I started using your function,
then I noticed it was very slow to execute.
Hence, I decided to show the process with a toy (linear) model.
The process remains the same.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize


def loss(parameters, model, x_ref, y_ref):
    # unpack multiple parameters
    m, q = parameters
    # compute prediction
    y_predicted = np.array([model(x, m, q) for x in x_ref])
    # compute error and use it as loss
    mse = ((y_ref - y_predicted) ** 2).mean()
    return mse


# load a dataset to fit a model
x_data = np.array([0.92, 0.925, 0.9535, 0.975, 0.9789, 1, 1.02, 1.043, 1.06, 1.078, 1.09])
y_data = np.array([150, 170, 269, 360, 377, 500, 583, 690, 761, 827, 840])
# normalise the data - input is already normalised
y_data_scaled = y_data / 1000

# create a model (linear, as an example) using handcrafted parameters, ex:(1,1)
linear_fun = lambda x, m, q: m * x + q
y_predicted = np.array([linear_fun(x, m=1, q=1) for x in x_data])

# create a function that given a model (linear_fun), a dataset(x,y) and the parameters, compute the error
parameters = (1, 1)
mse0 = loss(parameters, linear_fun, x_data, y_data_scaled)

# compute the parameters that minimize the loss (alias, the error between the data and the predictions of the model)
optimum = minimize(loss, x0=np.array([1, 1]), args=(linear_fun, x_data, y_data_scaled))

# compute the predictions with the optimized parameters
m = optimum['x'][0]
q = optimum['x'][1]
y_predicted_opt = np.array([linear_fun(x, m, q) for x in x_data])

# plot the raw data, the model with handcrafted guess and the model with optimized parameters
fig, ax = plt.subplots(1, 1)
ax.set_xlabel('input')
ax.set_ylabel('output predictions')
ax.plot(x_data, y_data_scaled, marker='o')
ax.plot(x_data, y_predicted, marker='*')
ax.plot(x_data, y_predicted_opt, marker='v')
ax.legend([
    "raw data points",
    "initial guess",
    "predictions with optimized parameters"
])

# Note1: good practise is to validate your model with a different set of data,
# respect to the one that you have used to find the parameters
# here, however, it is shown just the optimization procedure

# Note2: in your case you should use the HH_model instead of the linear_fun
# and I and Area_factor instead of m and q.

Output:

-- EDIT: To use the HH_model:

I went deeper in your code,
I tried few values for area and stimulus
and I executed a single run of HH_Model without taking the threshold.

Then, I checked the predicted dynamic of voltage (v):
the sequence is always diverging ( all values become nan after few steps )

if you have an initial guess for stimulus and area that could make the code to work, great.

if you have no idea of the order of magnitude of these parameters the unique solution I see is a grid search over them - just to find this initial guess.
however, it might take a very long time without guarantee of success.

given that the code is based on a physical model, I would suggest to:
1 - find pen and paper a reasonable values.
2 - check that this simulation works for these values.
3 - then, run the optimizer to find the minimum.

Or, worst case scenario, reverse engineer the code and find the value that makes the equation to converge

Here the refactored code:

import math
import numpy as np

# HH parameters
v_Rest = -65  # in mV
gNa = 1200  # in mS/cm^2
gK = 360  # in mS/cm^2
gL = 0.3 * 10  # in mS/cm^2
vNa = 115  # in mV
vK = -12  # in mV
vL = 10.6  # in mV
# Number of runs
c = 1  # in uF/cm^2


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


def predict_voltage(A, C, delay, dt, duration, stimulus, t_end):
    # 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.0005  # uA/(mS)^1/2

    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
        if delay / dt <= i <= (delay + duration) / dt:
            IStim = stimulus  # 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)[0]  # in ((uA / cm ^ 2) / (uF / cm ^ 2)) * ms == mV

        # stop simulation if it diverges
        if math.isnan(v[i + 1]):
            return [None]

    # adjust the voltage to the resting potential
    v = v + v_Rest
    return v


def HH_model(stimulus, area_factor, runs=1000):
    count = 0
    t_end = 10  # in ms
    delay = 0.1  # in ms
    duration = 0.1  # in ms
    dt = 0.0025  # in ms
    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 = (1 * 10 ** (-8)) * area_factor  # surface [cm^2]
    stimulus = stimulus
    C = c * A  # uF

    for j in range(0, runs):
        v = predict_voltage(A, C, delay, dt, duration, stimulus, t_end)
        if max(v[:] - v_Rest) > 60:
            count += 1
    return count

And the attempt to run one simulation:

import time
from ex_21.equations import c, predict_voltage

area_factor = 0.1
stimulus = 70

# input signal
count = 0
t_end = 10  # in ms
delay = 0.1  # in ms
duration = 0.1  # in ms
dt = 0.0025  # in ms

# geometry
d = 2  # diameter in um
r = d / 2  # Radius in um
l = 10  # Length of the compartment in  um
A = (1 * 10 ** (-8)) * area_factor  # surface [cm^2]
C = c * A  # uF

start = time.time()
voltage_dynamic = predict_voltage(A, C, delay, dt, duration, stimulus, t_end)
elapse = time.time() - start
print(voltage_dynamic)

Output:

[None]

Grid Search over function

1 Answers