Using full AP/CaT measures to fit the remaining conductances.
In [1]:
# PyABC imports
from pyabc import (ABCSMC, Distribution, RV,
History, MedianEpsilon)
from pyabc.populationstrategy import AdaptivePopulationSize
from pyabc.epsilon import MedianEpsilon
from pyabc.sampler import MulticoreEvalParallelSampler
In [2]:
# Custom imports
from ionchannelABC import (ion_channel_sum_stats_calculator,
IonChannelAcceptor,
IonChannelDistance,
EfficientMultivariateNormalTransition,
calculate_parameter_sensitivity,
plot_parameter_sensitivity,
plot_regression_fit,
plot_parameters_kde,
plot_distance_weights)
In [3]:
# Other necessary imports
import numpy as np
import subprocess
import pandas as pd
import io
import os
import tempfile
In [4]:
# Plotting imports
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%config Inline.Backend.figure_format = 'retina'
In [5]:
history = {'icat': History('sqlite:////storage/hhecm/cellrotor/chouston/results/icat-generic/hl-1_icat-generic.db'),
'ical': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ical-generic/hl-1_ical-generic.db'),
'ina': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ina-generic/hl-1_ina-generic3000.db'),
'ikr': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ikr-generic/hl-1_ikr-generic.db'),
'ikur': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ikur-generic/hl-1_ikur-generic.db'),
'ito': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ito-generic/hl-1_ito-generic.db'),
'iha': History('sqlite:////storage/hhecm/cellrotor/chouston/results/iha-generic/hl-1_iha-generic.db'),
'ik1': History('sqlite:////storage/hhecm/cellrotor/chouston/results/ik1-generic/hl-1_ik1-generic.db')}
In [6]:
run_ids = {'icat': 1, 'ical': 4, 'ina': 1, 'ikr': 2, 'ikur': 1, 'ito': 1, 'iha': 4, 'ik1': 1}
In [7]:
# Set to correct result (last ABC run in database).
for k,h in history.items():
h.id = run_ids[k]
In [8]:
# Create a prior from each channel's ABC posterior
from pyabc.transition import MultivariateNormalTransition
class rv_normal_fit(MultivariateNormalTransition):
def __init__(self, key):
self.key = key
super().__init__()
def fit(self, X, w):
X.columns = [self.key]
super().fit(pd.DataFrame(X), w)
def rvs(self):
sample = super().rvs()
return sample.values[0]
def pdf(self, x):
# wrap in pd.Series
x = pd.Series({self.key: x})
val = super().pdf(x)
return val
In [9]:
# Add to priors
prior = {}
for k, h in history.items():
df, w = h.get_distribution(m=0)
for key in df.keys():
dist = rv_normal_fit(k+'.'+key)
dist.fit(pd.DataFrame(df[key]), w)
prior[k+'.'+key] = dist
In [10]:
len(prior)
Out[10]:
In [11]:
from channels.hl1 import hl1 as model
In [12]:
#model.sample({})
In [13]:
# Targets from Dias, 2014 from HL1-6 myocyte
measurements = model.get_experiment_data()
obs = measurements.to_dict()['y']
exp = measurements.to_dict()['exp']
errs = measurements.to_dict()['errs']
In [14]:
measurements
Out[14]:
In [15]:
# Treat each entry as a separate experiment
for k, _ in exp.items():
exp[k] = k
In [16]:
# Add additional parameters to be fit using full cell
# measurements as uniform priors.
limits = {'ik1.g_K1': (0, 0.1),
'incx.k_NCX': (0, 20),
'icab.g_Cab': (0, 0.001),
'inab.g_Nab': (0, 0.01),
'inak.i_NaK_max': (0, 10)}
#'serca.V_max': (0, 10),
#'ca_diffusion.tau_tr': (0, 1000)}
#'ryanodine_receptors.k_RyR': (0, 0.1)}
for key, (a,b) in limits.items():
prior[key] = RV("uniform", a, b-a)
In [17]:
len(prior)
Out[17]:
In [18]:
prior = Distribution(**prior)
In [19]:
parameters = list(prior.keys())
In [20]:
distance_fn=IonChannelDistance(
obs=obs,
exp_map=exp,
err_bars=errs,
err_th=0.01)
In [21]:
sns.set_context('talk')
g = plot_distance_weights(model, distance_fn)
In [22]:
print(distance_fn.w)
In [23]:
g.savefig('results/extra-generic/dist_weights.pdf')
In [24]:
fitted, regression_fit, r2 = calculate_parameter_sensitivity(
model,
parameters,
distance_fn,
sigma=0.025,
n_samples=5000)
In [25]:
# Finding insensitive parameters
cutoff = 0.25
fitted_pivot = fitted.pivot(index='param',columns='exp')
insensitive_params = fitted_pivot[(abs(fitted_pivot['beta'][0])<cutoff) & (abs(fitted_pivot['beta'][1])<cutoff) &
(abs(fitted_pivot['beta'][2])<cutoff) & (abs(fitted_pivot['beta'][3])<cutoff) &
(abs(fitted_pivot['beta'][4])<cutoff) & (abs(fitted_pivot['beta'][5])<cutoff)].index.values
In [26]:
# SAVE RESULTS
import pickle
with open('insensitive.pkl','wb') as f:
pickle.dump(insensitive_params,f)
In [21]:
# LOAD RESULTS
import pickle
with open('insensitive.pkl','rb') as f:
insensitive_params = pickle.load(f)
In [27]:
len(insensitive_params)
Out[27]:
In [28]:
insensitive_params
Out[28]:
In [29]:
insensitive_params = np.array([p for p in insensitive_params if p not in limits.keys()])
In [30]:
len(insensitive_params)
Out[30]:
In [31]:
sns.set_context('talk')
grid = plot_parameter_sensitivity(fitted, plot_cutoff=0.2)
In [32]:
grid2 = plot_regression_fit(regression_fit, r2)
In [33]:
grid.savefig('results/extra-generic/sensitivity.pdf')
grid2.savefig('results/extra-generic/sensitivity_fit.pdf')
In [34]:
# Generate a sample from history items insensitive to full cell model measurements
def generate_sample(history, insensitive_params, n):
samples = [dict() for i in range(n)]
for k, h in history.items():
dist = h.get_distribution(m=0)
weights = dist[1]
ch_samples = dist[0] \
.sample(n, weights=weights, replace=True) \
.to_dict(orient='records')
ch_samples = [{k+'.'+key: value for key, value in ch.items()
if (k+'.'+key) in insensitive_params}
for ch in ch_samples]
for c, s in zip(ch_samples, samples):
s.update(c)
return samples
In [35]:
# Generate samples
n_samples = 10000
param_samples = generate_sample(history, insensitive_params, n_samples)
In [36]:
model.add_external_par_samples(param_samples)
In [37]:
# New prior without insensitive parameters
prior = {}
for k, h in history.items():
df, w = h.get_distribution(m=0)
for key in df.keys():
if (k+'.'+key) not in insensitive_params:
dist = rv_normal_fit(k+'.'+key)
dist.fit(pd.DataFrame(df[key]), w)
prior[k+'.'+key] = dist
In [38]:
limits = {'ik1.g_K1': (0, 0.1),
'incx.k_NCX': (0, 20),
'icab.g_Cab': (0, 0.001),
'inab.g_Nab': (0, 0.01),
'inak.i_NaK_max': (0, 10)}
#'serca.V_max': (0, 10),
#'ca_diffusion.tau_tr': (0, 1000)}
#'ryanodine_receptors.k_RyR': (0, 0.1)}
for key, (a,b) in limits.items():
prior[key] = RV("uniform", a, b-a)
In [39]:
prior = Distribution(**prior)
In [40]:
prior
Out[40]:
In [41]:
db_path = ('sqlite:///' +
os.path.join(tempfile.gettempdir(), "hl-1_extra-generic.db"))
print(db_path)
In [42]:
# Let's log all the sh!t
import logging
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
cv_logger = logging.getLogger('CV Estimation')
cv_logger.setLevel(logging.DEBUG)
In [43]:
from pyabc.sampler import SingleCoreSampler
from pyabc.populationstrategy import ConstantPopulationSize
In [44]:
abc = ABCSMC(models=model,
parameter_priors=prior,
distance_function=IonChannelDistance(
obs=obs,
exp_map=exp,
err_bars=errs,
err_th=0.01),
population_size=ConstantPopulationSize(5000),
#population_size=AdaptivePopulationSize(
# start_nr_particles=5000,
# mean_cv=0.2,
# max_population_size=5000,
# min_population_size=100),
summary_statistics=ion_channel_sum_stats_calculator,
transitions=EfficientMultivariateNormalTransition(scaling=1.0),
eps=MedianEpsilon(),#median_multiplier=1.5),
sampler=MulticoreEvalParallelSampler(n_procs=12),
#sampler=SingleCoreSampler(),
acceptor=IonChannelAcceptor())
In [45]:
abc_id = abc.new(db_path, obs)
In [ ]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=10, min_acceptance_rate=0.001)
In [ ]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=5, min_acceptance_rate=0.001)
In [ ]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=5, min_acceptance_rate=0.001)
In [ ]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=5, min_acceptance_rate=0.001)
In [ ]:
history = abc.run(minimum_epsilon=0.01, max_nr_populations=5, min_acceptance_rate=0.001)
In [48]:
#db_path = 'sqlite:///results/extra/hl-1_extra.db'
db_path = 'sqlite:////scratch/cph211/tmp/hl-1_extra-generic.db'
history = History(db_path)
history.all_runs()
Out[48]:
In [64]:
#history.id = 15 # run with err=0.1 and scaling=0.25
#history.id = 17 # run with err=0.01 and scaling=0.5
#history.id = 23 # run with new icat channel fitting using v_offset and no ion concentrations
#history.id = 24 # run with new icat channel fitting using v_offset, no ion concentrations and no APA
history.id = 1
In [65]:
evolution = history.get_all_populations()
grid = sns.relplot(x='t', y='epsilon', size='samples', data=evolution[evolution.t>=0])
#grid.savefig('results/extra/eps_evolution.pdf')
In [66]:
df_prior, w_prior = history.get_distribution(m=0,t=0)
In [67]:
df_extra = df_prior.loc[:,[k for k in prior.keys()]]
In [68]:
for c in df_extra.columns:
if c not in limits.keys():
limits[c] = (df_extra[c].min(), df_extra[c].max())
In [69]:
g = plot_parameters_kde(df_extra, w_prior, limits, aspect=12, height=0.6)
In [138]:
g.savefig('results/extra/parameters_kde_prior.pdf')
In [70]:
df, w = history.get_distribution(m=0)
In [71]:
df_extra = df.loc[:,[k for k in prior.keys()]]
In [72]:
limits
Out[72]:
In [73]:
g = plot_parameters_kde(df_extra, w, limits, aspect=12, height=0.6)
In [143]:
g.savefig('results/extra/parameters_kde_post.pdf')
In [74]:
# Generate parameter samples
n_samples = 100
df, w = history.get_distribution(m=0)
th_samples = df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records')
In [75]:
# Generate sim results samples
samples = pd.DataFrame({})
for i, th in enumerate(th_samples):
output = model.sample(pars=th)
output['sample'] = i
output['distribution'] = 'post'
samples = samples.append(output, ignore_index=True)
In [76]:
import scipy.stats as st
def evaluate_samples(samples):
for measurement in samples['x'].unique():
filtered = samples[samples.x==measurement]['y'].tolist()
print("{}, mean: {}".format(measurement, np.mean(filtered)))
print("{}, STD: {}".format(measurement, np.std(filtered)))
rv = st.rv_discrete(values=(filtered, [1/len(filtered),]*len(filtered)))
print("{}, median: {}".format(measurement, rv.median()))
print("{}, 95% CI: {}".format(measurement, rv.interval(0.95)))
In [77]:
evaluate_samples(samples)
In [78]:
from ionchannelABC import (Experiment, ExperimentStimProtocol)
In [79]:
# Create new experiment to generate traces
pulses = range(5, 45, 5)
stim_times = [1000000] + [2, 4998]*len(pulses)
stim_levels = [0]
for level in pulses:
stim_levels.append(level)
stim_levels.append(0)
#pace_time = np.linspace(0, 10000, 10000)
pace_time = np.linspace(0, sum(stim_times[1:]), sum(stim_times[1:]))
def trace_measurements(data_list):
import numpy as np
data = data_list[0]
if len(data_list) > 1:
for log in data_list[1:]:
data = data.extend(log)
simtime = data['environment.time']
simtime_min = min(simtime)
simtime = [t - simtime_min for t in simtime]
data_output = dict()
data['environment.time'] = simtime
for var in data:
data_output[var] = np.interp(pace_time, simtime, data[var])
return data_output
def unwrap(data, ind_var):
return data[0], True
trace_prot = ExperimentStimProtocol(stim_times, stim_levels,
#measure_index=range(len(stim_times)-4, len(stim_times)),
measure_index=range(1, len(stim_times)),
measure_fn=trace_measurements,
post_fn=unwrap)
dias_conditions = dict(T=305, Ca_o=1800, Na_o=1.4e5, K_o=4e3)
trace_exper = Experiment(trace_prot, None, dias_conditions)
In [80]:
from channels.hl1 import hl1 as model_testing
In [81]:
model_testing.add_experiments([trace_exper])
In [82]:
# Generate traces
traces = pd.DataFrame({})
for i, th in enumerate(th_samples[0:10]):
output = model_testing.sample(pars=th,exp_num=1)
output['sample'] = i
traces = traces.append(output, ignore_index=True)
In [145]:
import pickle
with open('traces.pkl','wb') as f:
pickle.dump(traces, f)
In [37]:
import pickle
with open('traces.pkl','rb') as f:
traces = pickle.load(f)
In [83]:
traces.x.unique()
Out[83]:
In [84]:
time = traces[traces.x=='environment.time'].y.tolist()[0]
In [85]:
#plotlims = (4950, 5400)
plotlims = [(4950+5000*i,5400+5000*i) for i in range(7)]
In [86]:
#idxlims = np.intersect1d(np.where(time >= plotlims[0]), np.where(time <= plotlims[1]))
idxlims = [np.intersect1d(np.where(time >= plotlims[i][0]), np.where(time <= plotlims[i][1]))
for i in range(7)]
In [87]:
#time_limited = time[idxlims[0]:idxlims[-1]]
time_limited = time[idxlims[0][0]:idxlims[0][-1]]
In [88]:
traces = traces.set_index(['x','sample'])
In [89]:
del traces['exp']
In [90]:
traces.reset_index(inplace=True)
In [91]:
rows = []
In [92]:
#_ = traces.apply(lambda row: [rows.append([row['x'], row['sample'], yi, ti])
# for (ti, yi) in zip(time_limited, row.y[idxlims[0]:idxlims[-1]])], axis=1)
for i in range(7):
_ = traces.apply(lambda row: [rows.append([row['x'], row['sample']+10*i, yi, ti])
for (ti, yi) in zip(time_limited, row.y[idxlims[i][0]:idxlims[i][-1]])], axis=1)
In [93]:
traces_new = pd.DataFrame(rows, columns=['x', 'sample', 'y', 'time']).set_index(['x', 'sample'])
In [94]:
traces_new.reset_index(inplace=True)
In [95]:
no_ap_samples = []
In [96]:
for i in range(max(traces_new['sample'])):
if max(traces_new[(traces_new['sample']==i) & (traces_new['x']=='membrane.V')]['y']) < 0:
no_ap_samples.append(i)
In [97]:
no_ap_samples
Out[97]:
In [98]:
for val in no_ap_samples:
traces_new = traces_new[~(traces_new['sample']==val)]
In [100]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='membrane.V')],
estimator=np.mean, ci='sd')
In [108]:
with sns.color_palette('gray'):
plot = sns.lineplot(x='time',y='y',hue='x',units='sample',
data=traces_new[(traces_new.x=='ikr.x')],
estimator=None,lw=0.5,alpha=0.5)
In [102]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='ca_conc.Ca_i')],
estimator=np.mean, ci='sd')
In [103]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='ina.i_Na')],
estimator=np.mean, ci='sd')
In [104]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='ical.i_CaL') | (traces_new.x=='icat.i_CaT')],
estimator=np.mean, ci='sd')
In [105]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='ikr.i_Kr') | (traces_new.x=='ikur.i_Kur') |
(traces_new.x=='ito.i_to')],
estimator=np.mean, ci='sd')
In [106]:
plot = sns.lineplot(x='time',y='y',hue='x',
data=traces_new[(traces_new.x=='ik1.i_K1') | (traces_new.x=='incx.i_NCX') |
(traces_new.x=='inak.i_NaK') | (traces_new.x=='iha.i_ha') |
(traces_new.x=='inab.i_Nab') | (traces_new.x=='icab.i_Cab')],
estimator=np.mean, ci='sd')
In [193]:
fig = plot.get_figure()
fig.savefig('results/extra/small_currents.pdf')
In [225]:
grid = sns.relplot(x='time', y='y', data=traces_new, row='x',
estimator=np.median, ci=95, kind='line',
facet_kws={'sharex': 'row',
'sharey': False})
In [49]:
samp_pivot = samples.pivot(index='sample', columns='x', values='y')
In [50]:
from sklearn.cluster import AffinityPropagation
In [51]:
af = AffinityPropagation()
In [52]:
af.fit(samp_pivot.values)
Out[52]:
In [53]:
centers = pd.DataFrame(af.cluster_centers_, columns=samp_pivot.columns)
In [55]:
af.cluster_centers_indices_
Out[55]:
In [56]:
sns.distplot(af.labels_, bins=len(af.cluster_centers_))
Out[56]: