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)
In [7]:
'VClamp' in sampleData, 'IClamp' in sampleData
Out[7]:
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()