Parameter Extraction for EGL-19 Ion Channel


In [1]:
"""
Example of using cwFitter to generate a HH model for EGL-19 Ca2+ ion channel
Based on experimental data from doi:10.1083/jcb.200203055
"""
import os.path
import sys
sys.path.append('..')
sys.path.append('../..')
sys.path.append('../../..')

import django

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
from neurotune import optimizers

from channelworm.fitter import *

# Setup access to the Django database
os.environ.setdefault(
    "DJANGO_SETTINGS_MODULE",
    "channelworm.web_app.settings"
)
django.setup()

from channelworm.ion_channel.models import GraphData

In [2]:
def IV_act(V,g,Vhalf,k,a_power,e_rev):
    return g * (1/(1 + np.exp((Vhalf - V)/k)))**int(a_power) * (V - e_rev)

# if __name__ == '__main__':

userData = dict()

cwd=os.getcwd()
csv_path = os.path.dirname(cwd)+'/examples/egl-19-data/egl-19-IClamp-IV.csv'
ref = {'fig':'2B','doi':'10.1083/jcb.200203055'}
x_var = {'type':'Voltage','unit':'V','toSI':1}
y_var = {'type':'Current','unit':'A/F','toSI':1}
IV = {'ref':ref,'csv_path':csv_path,'x_var':x_var,'y_var':y_var}
userData['samples'] = {'IV':IV}

In [3]:
myInitiator = initiators.Initiator(userData)

In [4]:
myInitiator = initiators.Initiator(userData)
sampleData = myInitiator.get_sample_params()
bio_params = myInitiator.get_bio_params()
sim_params = myInitiator.get_sim_params(type='IClamp')
# sim_params = myInitiator.get_sim_params()
sim_params = sim_params['IC']
sim_params['ca_con'] = 1e-6

myEvaluator = evaluators.Evaluator(sampleData,sim_params,bio_params)

"""
Using a population size 20 times larger than the number of free parameters 
is what is suggested by Gurkiewicz and Korngreen.  
"""

candidates = optimizers.CustomOptimizerA(bio_params['max_val_channel'],
                                             bio_params['min_val_channel'],
                                             myEvaluator,
                                             population_size=300, #20 times larger than free parameters
                                             max_evaluations=600,
                                             num_selected=2,
                                             num_offspring=15,
                                             num_elites=1,
                                             mutation_rate=0.05,
                                             maximize = False,
                                             seeds=None,
                                             verbose=True)

In [5]:
# best_candidate = candidates.optimize(do_plot=True, seed=1234)
#"""
best_candidate = [189.66504757196532, 0.04619592278775828, -0.0015715032129984834, 0.03232782689368996,
   0.0009038799935426481, -0.0006996189007248855, 0.0002076054785033701, 0.5361776032113692,
   2.0, 1.0, 2.9942088447494227e-07, 0.18712673917281542, -1.1396759086697654e-07,
   0.014145060464901655, 1.0]
#"""
best_candidate_params = dict(zip(bio_params['channel_params'],best_candidate))
cell_var = dict(zip(bio_params['cell_params'],bio_params['val_cell_params']))
mySimulator = simulators.Simulator(sim_params,best_candidate_params,cell_var,bio_params['gate_params'])
bestSim = dict()

In [6]:
# The I/V plot could come from either VClamp or IClamp (VClamp is preferred as is more common)

if ('VClamp' in sampleData) or (('IV' in sampleData) and (('VClamp' and 'IClamp') not in sampleData)):
  bestSim.update({'VClamp':{}})
  bestSim['VClamp']['t'],bestSim['VClamp']['I'],bestSim['VClamp']['V_max'],bestSim['VClamp']['I_max'] = mySimulator.VClamp()
if 'IClamp' in sampleData:
  bestSim.update({'IClamp':{}})
  bestSim['IClamp']['t'],bestSim['IClamp']['V'],bestSim['IClamp']['V_max'],bestSim['IClamp']['I_max'] = mySimulator.IClamp()

myModelator = modelators.Modelator(bio_params,sim_params).compare_plots(sampleData,bestSim,show=True)


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-6-b02b170e77de> in <module>()
      3 if ('VClamp' in sampleData) or (('IV' in sampleData) and (('VClamp' and 'IClamp') not in sampleData)):
      4   bestSim.update({'VClamp':{}})
----> 5   bestSim['VClamp']['t'],bestSim['VClamp']['I'],bestSim['VClamp']['V_max'],bestSim['VClamp']['I_max'] = mySimulator.VClamp()
      6 if 'IClamp' in sampleData:
      7   bestSim.update({'IClamp':{}})

AttributeError: 'Simulator' object has no attribute 'VClamp'

In [7]:
'VClamp' in sampleData, 'IClamp' in sampleData


Out[7]:
(False, False)

In [ ]:
# Fitting to the I/V curve and optimizing parameters
# According to the literature, the I/V plot coming from steady state currents
# Only activation expressions will be considered
vData = np.arange(-0.040, 0.080, 0.001)
p0 = [best_candidate_params['g'],best_candidate_params['v_half_a'],best_candidate_params['k_a'],best_candidate_params['a_power'],best_candidate_params['e_rev']]

Vsample = np.asarray(sampleData['IV']['V'])
Isample = np.asarray(sampleData['IV']['I'])
popt,pcov = curve_fit(IV_act, Vsample,Isample,p0)
Iopt = IV_act(vData,popt[0],popt[1],popt[2],popt[3],popt[4])


print(p0)
print(popt)

In [ ]:
if 'VClamp' in bestSim:
  model_plot, = plt.plot([x*1e3 for x in bestSim['VClamp']['V_max']],bestSim['VClamp']['I_max'], label = 'simulated using GA')
else:
  model_plot, = plt.plot([x*1e3 for x in bestSim['IClamp']['V_max']],bestSim['IClamp']['I_max'])
  sample, = plt.plot([x*1e3 for x in sampleData['IV']['V']],sampleData['IV']['I'], 'y', label = 'sample data')
  # sim, =plt.plot(vData,I, label = 'simulated_curve')
  opt, =plt.plot([x*1e3 for x in vData],Iopt, 'r', label = 'optimized with GA and  minimization function')
  plt.legend([model_plot,sample,opt])
  # plt.legend([model_plot,sample,opt,sim])
  # plt.legend([sample,opt])
  plt.title("The Best Model fitted to data using GA and  minimization function")
  plt.ylabel('I (A/F)')
  plt.xlabel('V (mV)')
  plt.show()