Analysis of results from rapidly-activating rectifier Potassium channel ABC calibration

This notebook compares the behaviour of the $I_\text{Kur}$ current in Nygren, Courtemanche and standardised formulations. Nygren terminology $I_\text{sus}$ is used throughout the notebook.



In [1]:

    
import pandas as pd
import numpy as np
import scipy.stats as stats

import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt

from pyabc import History
from pyabc.weighted_statistics import weighted_mean, weighted_median

from ionchannelABC.visualization import plot_variables
from ionchannelABC.utils import weighted_cv









    



INFO:myokit:Loading Myokit version 1.28.3

Load pyABC results databases for each run

Assumes that the results databases are stored in a subdirectory called `results'.

Loads history files and converts to dataframes and weights (see pyABC documentation).

Nygren



In [2]:

    
h_nyg_r_original = History('sqlite:///results/nygren/isus/original/nygren_isus_rgate_original.db')
h_nyg_s_original = History('sqlite:///results/nygren/isus/original/nygren_isus_sgate_original.db')



In [3]:

    
h_nyg_r_unified = History('sqlite:///results/nygren/isus/unified/nygren_isus_rgate_unified.db')



In [4]:

    
df_nyg_r_original, w_nyg_r_original = h_nyg_r_original.get_distribution(m=0)
df_nyg_s_original, w_nyg_s_original = h_nyg_s_original.get_distribution(m=0)



In [5]:

    
df_nyg_r_unified, w_nyg_r_unified = h_nyg_r_unified.get_distribution(m=0)

Courtemanche



In [6]:

    
h_cou_a_original = History('sqlite:///results/courtemanche/isus/original/courtemanche_isus_agate_original.db')
h_cou_i_original = History('sqlite:///results/courtemanche/isus/original/courtemanche_isus_igate_original.db')



In [7]:

    
h_cou_i_unified = History('sqlite:///results/courtemanche/isus/unified/courtemanche_isus_igate_unified.db')



In [8]:

    
df_cou_a_original, w_cou_a_original = h_cou_a_original.get_distribution(m=0)
df_cou_i_original, w_cou_i_original = h_cou_i_original.get_distribution(m=0)



In [9]:

    
df_cou_i_unified, w_cou_i_unified = h_cou_i_unified.get_distribution(m=0)

Standardised



In [10]:

    
h_sta = History('sqlite:///results/standardised/isus/standardised_isus.db')



In [11]:

    
df_sta, w_sta = h_sta.get_distribution(m=0)

Compare number of gating parameters in models

Generate bar plot using manually input numbers. See models and publications for source.## Compare number of parameters



In [12]:

    
data = [['N',6,7],['C',13,7],['S',4,4]]



In [13]:

    
df = pd.DataFrame(data, columns = ['model', 'act', 'inact'])



In [14]:

    
df.inact = df.act+df.inact



In [15]:

    
sns.set(style="ticks")
sns.set_context('talk')
sns.set_color_codes("pastel")

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plt.subplots(figsize=(6,3))

sns.barplot(x="inact", y="model", hue="model", data=df,
            palette="pastel", dodge=False)
            #label="inact", color="b")

sns.set_color_codes("muted")
sns.barplot(x="act", y="model", hue="model", data=df,
            palette="muted", dodge=False)
            #label="act", color="r")

handles, labels = ax.get_legend_handles_labels()
order = []
ax.legend([handles[idx] for idx in order],[labels[idx] for idx in order],
          ncol=2, loc="lower right", frameon=False)
ax.set(xlabel="number of gating parameters")
ax.xaxis.grid(True)
ax.set(ylabel="")
sns.despine(trim=True, left=True)



In [16]:

    
#f.savefig('figures/isus/param_num_compare.pdf')

Compare calibration between original and unified datasets

Plot relative standard deviation (RSD) of parameter posteriors and gating functions for Nygren and Courtemanche using the original dataset and unified datasets.

Note: we use the term RSD instead of CV (coefficient of variation) to avoid confusion with conduction velocity.

Calculate RSD values for Nygren and Courtemanche models



In [17]:

    
rsd_nyg_r_original = weighted_cv(df_nyg_r_original, w_nyg_r_original, sample_size=len(df_nyg_r_original))
rsd_nyg_s_original = weighted_cv(df_nyg_s_original, w_nyg_s_original, sample_size=len(df_nyg_s_original))

rsd_nyg_original = pd.concat([rsd_nyg_r_original, rsd_nyg_s_original])
model = ['N',]*len(rsd_nyg_original)
data = ['ORIGINAL',]*len(rsd_nyg_original)
frame = {'model': model, 'data': data, 'RSD': rsd_nyg_original}
nyg_original = pd.DataFrame(frame)



In [18]:

    
rsd_nyg_r_unified = weighted_cv(df_nyg_r_unified, w_nyg_r_unified, sample_size=len(df_nyg_r_unified))

rsd_nyg_unified = pd.concat([rsd_nyg_r_unified, rsd_nyg_s_original])
model = ['N',]*len(rsd_nyg_unified)
data = ['UNIFIED',]*len(rsd_nyg_unified)
frame = {'model': model, 'data': data, 'RSD': rsd_nyg_unified}
nyg_unified = pd.DataFrame(frame)



In [19]:

    
rsd_cou_a_original = weighted_cv(df_cou_a_original, w_cou_a_original, sample_size=len(df_cou_a_original))
rsd_cou_i_original = weighted_cv(df_cou_i_original, w_cou_i_original, sample_size=len(df_cou_i_original))

rsd_cou_original = pd.concat([rsd_cou_a_original, rsd_cou_i_original])
model = ['C',]*len(rsd_cou_original)
data = ['ORIGINAL',]*len(rsd_cou_original)
frame = {'model': model, 'data': data, 'RSD': rsd_cou_original}
cou_original = pd.DataFrame(frame)



In [20]:

    
rsd_cou_i_unified = weighted_cv(df_cou_i_unified, w_cou_i_unified, sample_size=len(df_cou_i_unified))

rsd_cou_unified = pd.concat([rsd_cou_a_original, rsd_cou_i_unified])
model = ['C',]*len(rsd_cou_unified)
data = ['UNIFIED',]*len(rsd_cou_unified)
frame = {'model': model, 'data': data, 'RSD': rsd_cou_unified}
cou_unified = pd.DataFrame(frame)



In [21]:

    
rsd_compare = pd.concat([nyg_original, nyg_unified, cou_original, cou_unified], sort=False)



In [22]:

    
sns.set(style="ticks")
sns.set_context("poster")

# Initialize the figure
f, ax = plt.subplots(figsize=(4, 6))

# Plot the boxplot summary of RSD
sns.boxplot(x="model", y="RSD", hue="data", data=rsd_compare,
            palette="pastel", whis="range")

# Add in points to show each observation
sns.swarmplot(x="model", y="RSD", hue="data", data=rsd_compare,
              linewidth=0.5, size=8, dodge=True)

# Tweak the visual presentation
ax.yaxis.grid(True)
ax.set(ylabel="RSD")
ax.set_yscale('log')
ax.set(xlabel="")
sns.despine(trim=True, bottom=True)
ax.get_legend().remove()



In [23]:

    
#f.savefig('figures/isus/rsd_compare.pdf')

Compare underlying gating functions for each model

Sample from posterior particle population and plot resulting gating functions



In [24]:

    
N = 100



In [25]:

    
nyg_par_samples_r_original = df_nyg_r_original.sample(n=N, weights=w_nyg_r_original, replace=True)
nyg_par_samples_r_original = nyg_par_samples_r_original.set_index([pd.Index(range(N))])
nyg_par_samples_s_original = df_nyg_s_original.sample(n=N, weights=w_nyg_s_original, replace=True)
nyg_par_samples_s_original = nyg_par_samples_s_original.set_index([pd.Index(range(N))])

nyg_par_samples_original = (pd.concat([nyg_par_samples_r_original, 
                                       nyg_par_samples_s_original],axis=1)
                            .to_dict(orient='records'))



In [26]:

    
nyg_par_samples_r_unified = df_nyg_r_unified.sample(n=N, weights=w_nyg_r_unified, replace=True)
nyg_par_samples_r_unified = nyg_par_samples_r_unified.set_index([pd.Index(range(N))])

nyg_par_samples_unified = (pd.concat([nyg_par_samples_r_unified, 
                                      nyg_par_samples_s_original],axis=1)
                           .to_dict(orient='records'))



In [27]:

    
v = np.arange(-100, 50, 0.5)

nyg_par_map = {'ri': 'isus.r_inf',
               'rt': 'isus.tau_r',
               'si': 'isus.s_inf',
               'st': 'isus.tau_s'}



In [28]:

    
from ionchannelABC.visualization import plot_variables
sns.set_context('poster')
sns.set_palette('deep')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plot_variables(v, [nyg_par_map, nyg_par_map], 
                       ['models/nygren_isus.mmt', 'models/nygren_isus_.mmt'], 
                       [nyg_par_samples_original, nyg_par_samples_unified],
                       original=[True,False],
                       figshape=(2,2))

ax[0][0].set_ylabel(r'$a$')
ax[1][0].set_ylabel(r'$i$')

ax[0][0].set_title(r'$\bar{x}$')
ax[0][1].set_title(r'$\tau_x$ (ms)')

for a in ax[:,1].flatten():
    a.set_ylabel('')
    
for a in ax[1][:]:
    a.set_xlabel('voltage (mV)')

for a in ax[:,0]:
    a.set_ylim((-0.05, 1.05))
    
plt.tight_layout()









    



/scratch/cph211/miniconda3/envs/ionchannelABC/lib/python3.7/site-packages/pandas/core/frame.py:6692: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version
of pandas will change to not sort by default.

To accept the future behavior, pass 'sort=False'.

To retain the current behavior and silence the warning, pass 'sort=True'.

  sort=sort)



In [29]:

    
#f.savefig('figures/isus/nyg_gating_functions.pdf')



In [30]:

    
cou_par_samples_a_original = df_cou_a_original.sample(n=N, weights=w_cou_a_original, replace=True)
cou_par_samples_a_original = cou_par_samples_a_original.set_index([pd.Index(range(N))])
cou_par_samples_i_original = df_cou_i_original.sample(n=N, weights=w_cou_i_original, replace=True)
cou_par_samples_i_original = cou_par_samples_i_original.set_index([pd.Index(range(N))])

cou_par_samples_original = pd.concat([cou_par_samples_a_original, cou_par_samples_i_original],axis=1).to_dict(orient='records')



In [31]:

    
cou_par_samples_i = df_cou_i_unified.sample(n=N, weights=w_cou_i_unified, replace=True)
cou_par_samples_i = cou_par_samples_i.set_index([pd.Index(range(N))])
cou_par_samples_unified = pd.concat([cou_par_samples_a_original, cou_par_samples_i],axis=1).to_dict(orient='records')



In [32]:

    
cou_par_map = {'ri': 'isus.a_inf',
               'rt': 'isus.tau_a',
               'si': 'isus.i_inf',
               'st': 'isus.tau_i'}



In [33]:

    
from ionchannelABC.visualization import plot_variables
sns.set_context('poster')
sns.set_palette('deep')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plot_variables(v, [cou_par_map, cou_par_map], 
                       ['models/courtemanche_isus.mmt', 'models/courtemanche_isus_.mmt'], 
                       [cou_par_samples_original, cou_par_samples_unified],
                       original=[True,False],
                       figshape=(2,2))

ax[0][0].set_ylabel(r'$a$')
ax[1][0].set_ylabel(r'$i$')

ax[0][0].set_title(r'$\bar{x}$')
ax[0][1].set_title(r'$\tau_x$ (ms)')

for a in ax[:,1].flatten():
    a.set_ylabel('')
    
for a in ax[1][:]:
    a.set_xlabel('voltage (mV)')

for a in ax[:,0]:
    a.set_ylim((-0.05, 1.05))

ax[1,1].set_ylim((-0.05*25000, 1.05*25000))
    
plt.tight_layout()



In [34]:

    
#f.savefig('figures/isus/cou_gating_functions.pdf')

Compare all models calibrated to unified dataset

RSD of parameter posteriors



In [35]:

    
rsd_sta = weighted_cv(df_sta, w_sta, sample_size=len(df_sta))
model = ['S',]*len(rsd_sta)
frame = {'model': model, 'RSD': rsd_sta}
sta = pd.DataFrame(frame)



In [36]:

    
rsd_unified = pd.concat([nyg_unified, cou_unified, sta], sort=False)



In [37]:

    
sns.set(style="ticks")
sns.set_context('talk')

# Initialize the figure
f, ax = plt.subplots(figsize=(6, 3))

sns.boxplot(x="RSD", y="model", data=rsd_unified,
            palette="pastel", whis="range")

# Add in points to show each observation
sns.swarmplot(x="RSD", y="model", data=rsd_unified,
              linewidth=0.5, size=8, dodge=True)

# Tweak the visual presentation
ax.xaxis.grid(True)
ax.set(xlabel="RSD")
ax.set_xscale("log")
ax.set(ylabel="")
sns.despine(trim=True, left=True)



In [38]:

    
#f.savefig('figures/isus/rsd_compare_unified.pdf')



In [39]:

    
stats.mannwhitneyu(rsd_unified[rsd_unified.model=='N'].RSD, 
                   rsd_unified[rsd_unified.model=='S'].RSD)









    Out[39]:





MannwhitneyuResult(statistic=45.0, pvalue=0.31891519499830134)



In [40]:

    
stats.mannwhitneyu(rsd_unified[rsd_unified.model=='C'].RSD, 
                   rsd_unified[rsd_unified.model=='S'].RSD)









    Out[40]:





MannwhitneyuResult(statistic=32.0, pvalue=0.00785462752352008)



In [41]:

    
stats.mannwhitneyu(rsd_unified[rsd_unified.model=='N'].RSD, 
                   rsd_unified[rsd_unified.model=='C'].RSD)









    Out[41]:





MannwhitneyuResult(statistic=115.0, pvalue=0.2965894121668926)

Residuals to experimental data



In [42]:

    
N = 100



In [43]:

    
eps_nyg_r = h_nyg_r_unified.get_weighted_distances()
eps_nyg_s = h_nyg_s_original.get_weighted_distances()



In [44]:

    
eps_nyg = np.array([])
eps_nyg = (np.array(eps_nyg_r.sample(n=100,axis=0,weights=eps_nyg_r.w,replace=True).distance) +
           np.array(eps_nyg_s.sample(n=100,axis=0,weights=eps_nyg_s.w,replace=True).distance))
eps_nyg = pd.DataFrame({'model': 'N', 'eps': eps_nyg, 'exp': 'all'})



In [45]:

    
eps_cou_a = h_cou_a_original.get_weighted_distances()
eps_cou_i = h_cou_i_unified.get_weighted_distances()
eps_cou = np.array([])
eps_cou = (np.array(eps_cou_a.sample(n=100,axis=0,weights=eps_cou_a.w,replace=True).distance) +
           np.array(eps_cou_i.sample(n=100,axis=0,weights=eps_cou_i.w,replace=True).distance))
eps_cou = pd.DataFrame({'model': 'C', 'eps': eps_cou, 'exp': 'all'})



In [46]:

    
eps_sta = h_sta.get_weighted_distances()
eps_sta = eps_sta.sample(n=100,axis=0,weights=eps_sta.w,replace=True).distance
eps_sta = pd.DataFrame({'model': 'S', 'eps': eps_sta, 'exp': 'all'})



In [47]:

    
eps = pd.concat([eps_nyg,eps_cou,eps_sta])



In [48]:

    
eps['normalised'] = (eps['eps']-eps['eps'].min())/(eps['eps'].max()-eps['eps'].min())



In [49]:

    
sns.set_context("talk")
sns.set_style("whitegrid")

mpl.rcParams["font.size"] = 14
mpl.rcParams["legend.fontsize"] = 14

f, ax = plt.subplots(figsize=(2, 3))
sns.despine(bottom=True, left=True)

sns.stripplot(x='exp', y='normalised', hue='model', data=eps, 
              dodge=True, jitter=True, alpha=.25, zorder=1,
              palette='deep', ax=ax)

sns.pointplot(x="exp", y="normalised", hue="model",
              data=eps,
              estimator=np.median,
              dodge=.532, join=False, palette="dark",
              markers="d", scale=.75, ci=None)

ax.get_legend().remove()
ax.set_xlabel('experiment number')
ax.set_ylabel('normalised residual')









    Out[49]:





Text(0, 0.5, 'normalised residual')



In [50]:

    
#f.savefig('figures/isus/eps_all_stripplot.pdf')



In [51]:

    
sns.set_context('talk')
sns.set_style("ticks")

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plt.subplots(figsize=(6,3))

sns.barplot(x="eps", y="model", data=eps, palette='pastel')

ax.xaxis.grid(True)
ax.set(ylabel="")
ax.set(xlabel="converged residual")
sns.despine(trim=True, left=True)



In [52]:

    
#f.savefig('figures/isus/eps_all_barplot.pdf')

Below we sample from the posterior particle populations and run the calibration experiments to be able to calculate a per-experiment epsilon value.



In [53]:

    
N = 100

nyg_par_samples_r = df_nyg_r_unified.sample(n=N, weights=w_nyg_r_unified, replace=True)
nyg_par_samples_r = nyg_par_samples_r.set_index([pd.Index(range(N))])
nyg_par_samples_r = nyg_par_samples_r.to_dict(orient='records')
nyg_par_samples_s = df_nyg_s_original.sample(n=N, weights=w_nyg_s_original, replace=True)
nyg_par_samples_s = nyg_par_samples_s.set_index([pd.Index(range(N))])
nyg_par_samples_s = nyg_par_samples_s.to_dict(orient='records')

cou_par_samples_a = df_cou_a_original.sample(n=N, weights=w_cou_a_original, replace=True)
cou_par_samples_a = cou_par_samples_a.set_index([pd.Index(range(N))])
cou_par_samples_a = cou_par_samples_a.to_dict(orient='records')
cou_par_samples_i = df_cou_i_unified.sample(n=N, weights=w_cou_i_unified, replace=True)
cou_par_samples_i = cou_par_samples_i.set_index([pd.Index(range(N))])
cou_par_samples_i = cou_par_samples_i.to_dict(orient='records')

std_par_samples = df_sta.sample(n=N, weights=w_sta, replace=True)
std_par_samples = std_par_samples.set_index([pd.Index(range(N))])
std_par_samples = std_par_samples.to_dict(orient='records')



In [54]:

    
from ionchannelABC.experiment import setup
from ionchannelABC.distance import IonChannelDistance



In [58]:

    
from experiments.isus_wang import (wang_act,
                                   wang_act_kin,
                                   wang_inact)
from experiments.isus_courtemanche import (courtemanche_inact_kin,
                                           courtemanche_deact)
from experiments.isus_firek import (firek_inact)
from experiments.isus_nygren import (nygren_inact_kin,
                                     nygren_rec)

experiments = [wang_act,
               wang_act_kin,
               courtemanche_deact,
               firek_inact,
               nygren_inact_kin,
               nygren_rec]
models = ['models/nygren_isus_temp_adj.mmt',
          'models/courtemanche_isus_temp_adj.mmt',
          'models/standardised_isus.mmt']
par_samples = [[nyg_par_samples_r,nyg_par_samples_s], [cou_par_samples_a,cou_par_samples_i], std_par_samples]
names = ['N','C','S']



In [59]:

    
def experiment_dist(par_samples, modelfile, experiment):
    eps = []
    observ, model, sum_stats = setup(modelfile,
                                     experiment)
    obs = observ.to_dict()['y']
    obs = {str(k): v for k, v in obs.items()}
    dist = IonChannelDistance(exp_id=list(observ.exp_id),
                              variance=list(observ.variance),
                              delta=0.05)
    for sample in par_samples:
        eps.append(dist(sum_stats(model(sample)), obs, None))    
    return eps

The below cell calculates the per-experiment distance and may take some time to run depending on number of samples.



In [60]:

    
df = pd.DataFrame({})
mask = [[0, 0, 0, 1, 1, 1], [0, 0, 0, 1, 1, 1], None]
for j, exp in enumerate(experiments):
    df_ = pd.DataFrame({})
    for i, m in enumerate(models):
        if mask[i] is not None:
            eps = np.array(
                experiment_dist(par_samples[i][mask[i][j]],
                                m,
                                exp)
                )
        else:
            eps = np.array(
                experiment_dist(par_samples[i],
                                m,
                                exp)
                )
        df_ = df_.append(pd.DataFrame({'model': names[i],
                                       'exp': str(j),
                                       'eps': eps[~np.isinf(eps)]}))
    eps_max = df_['eps'].max()
    eps_min = df_['eps'].min()
    df_['eps'] = (df_['eps'] - eps_min)/(eps_max-eps_min)
    df = df.append(df_)



In [61]:

    
sns.set_context('talk')
sns.set_style('whitegrid')

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

fig, ax = plt.subplots(figsize=(12, 3))
sns.despine(bottom=True, left=True)

sns.stripplot(x='exp', y='eps', hue='model', data=df, 
            dodge=True, jitter=True, alpha=.25, zorder=1,
              palette='deep', ax=ax)

sns.pointplot(x="exp", y="eps", hue="model",
              data=df,
              estimator=np.median,
              dodge=.532, join=False, palette="dark",
              markers="d", scale=.75, ci=None)

handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[3:], labels[3:], title="model",
          handletextpad=0, columnspacing=1,
          loc="best", ncol=3, frameon=True)
ax.set_xlabel('experiment number')
ax.set_ylabel('normalised residual')









    Out[61]:





Text(0, 0.5, 'normalised residual')



In [62]:

    
#fig.savefig('figures/isus/eps_per_exp_stripplot.pdf')

Summary statistics compared to experimental data



In [63]:

    
from ionchannelABC.visualization import plot_sim_results



In [64]:

    
from experiments.isus_wang import (wang_act_and_kin,
                                   wang_inact)
from experiments.isus_courtemanche import (courtemanche_inact_kin,
                                           courtemanche_deact)
from experiments.isus_firek import (firek_inact)
from experiments.isus_nygren import (nygren_inact_kin,
                                     nygren_rec)



In [67]:

    
sns.set(style="ticks")
sns.set_context('talk')

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

g = plot_sim_results(['models/nygren_isus.mmt',
                      'models/nygren_isus.mmt',
                      'models/courtemanche_isus.mmt',
                      'models/courtemanche_isus.mmt',
                      'models/standardised_isus.mmt'],
                     wang_act_and_kin,
                    courtemanche_deact,
                    firek_inact,
                    nygren_inact_kin,
                    nygren_rec,
                     temp_match_model = 4,
                     masks=[[(0,1),2]+[None,]*3,
                            [None,]*2+[3,4,5],
                            [(0,1),2]+[None,]*3,
                            [None,]*2+[3,4,5],
                            None],
                     df=[df_nyg_r_unified,
                         df_nyg_s_original,
                         df_cou_a_original,
                         df_cou_i_unified,
                         df_sta], 
                     w=[w_nyg_r_unified,
                        w_nyg_s_original,
                        w_cou_a_original,
                        w_cou_i_unified,
                        w_sta])

ylabels = ["steady-state activation", "activation τ (ms)", "deactivation τ (ms)", "steady-state inactivation",
           "inactivation τ (ms)", "recovery τ (ms)"]
for ax in g.axes.flatten():
    ax.set_xlabel("voltage (mV)")
for ax, yl in zip(g.axes.flatten(), ylabels):
    ax.set_ylabel(yl)
for ax in g.axes.flatten():
    ax.set_title("")

plt.tight_layout()



In [72]:

    
#g.savefig('figures/isus/compare_summary_statistics.pdf')

Compare SS and tau values



In [73]:

    
nyg = 'models/nygren_isus_temp_adj.mmt'
cou ='models/courtemanche_isus_temp_adj.mmt'
std = 'models/standardised_isus.mmt'



In [74]:

    
v = np.arange(-80, 50, 0.5)



In [75]:

    
std_par_map = {'ri': 'isus.r_ss',
               'si': 'isus.s_ss',
               'rt': 'isus.tau_r',
               'st': 'isus.tau_s'}



In [76]:

    
N = 100

nyg_par_samples_r = df_nyg_r_unified.sample(n=N, weights=w_nyg_r_unified, replace=True)
nyg_par_samples_r = nyg_par_samples_r.set_index([pd.Index(range(N))])
nyg_par_samples_s = df_nyg_s_original.sample(n=N, weights=w_nyg_s_original, replace=True)
nyg_par_samples_s = nyg_par_samples_s.set_index([pd.Index(range(N))])
nyg_par_samples = (pd.concat([nyg_par_samples_r,
                              nyg_par_samples_s],axis=1).to_dict(orient='records'))

cou_par_samples_a = df_cou_a_original.sample(n=N, weights=w_cou_a_original, replace=True)
cou_par_samples_a = cou_par_samples_a.set_index([pd.Index(range(N))])
cou_par_samples_i = df_cou_i_unified.sample(n=N, weights=w_cou_i_unified, replace=True)
cou_par_samples_i = cou_par_samples_i.set_index([pd.Index(range(N))])
cou_par_samples = pd.concat([cou_par_samples_a,cou_par_samples_i],axis=1).to_dict(orient='records')

std_par_samples = df_sta.sample(n=N, weights=w_sta, replace=True)
std_par_samples = std_par_samples.set_index([pd.Index(range(N))])
std_par_samples = std_par_samples.to_dict(orient='records')



In [78]:

    
sns.set(style='ticks')
sns.set_context('poster')

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plot_variables(v, [nyg_par_map, cou_par_map, std_par_map], 
                       [nyg, cou, std], 
                       [nyg_par_samples, cou_par_samples, std_par_samples],
                       original=[True, True, False],
                       figshape=(2,2))

ax[0][0].set_title('steady state')
ax[0][1].set_title('time constant (ms)')
ax[0][0].set_ylabel('Activation')
ax[1][0].set_ylabel('Inactivation')

for a in ax[:,0]:
    a.set_ylim((-0.05, 1.05))

for a in ax[:,1]:
    a.set_ylabel('')
    
for a in ax[1,:]:
    a.set_xlabel('voltage (mV)')
    
ax[1][1].set_ylim([-0.05*100+200, 1.05*100+200])

    
for a in ax.flatten():
    x0,x1 = a.get_xlim()
    y0,y1 = a.get_ylim()
    a.set_aspect(abs(x1-x0)/abs(y1-y0))
    
plt.tight_layout()



In [79]:

    
#f.savefig('figures/isus/compare_gating_functions.pdf')

Compare final pulse after pulse train



In [80]:

    
import myokit



In [81]:

    
nsteps = 100
period = 10000
tstep = 1000
vhold = -50
vstep = -10

pulsetrain = myokit.pacing.steptrain([vstep,]*nsteps,
                                     vhold,
                                     period-tstep,
                                     tstep)



In [82]:

    
# nygren_ina_full has had tau values adjusted for temperature to 310K
models = [myokit.load_model(modelfile) for modelfile in ['models/nygren_isus_temp_adj.mmt',
                                                         'models/courtemanche_isus_temp_adj.mmt',
                                                         'models/standardised_isus.mmt']]
for m in models:
    pace = m.get('membrane.V')
    if pace.binding() != 'pace':
        if pace.is_state():
            pace.demote()
        pace.set_rhs(0)
        pace.set_binding('pace')



In [83]:

    
sims = [myokit.Simulation(m, pulsetrain) for m in models]



In [84]:

    
par_samples = [nyg_par_samples, cou_par_samples, std_par_samples]



In [85]:

    
import time

all_samples = pd.DataFrame({})
for i, s in enumerate(sims):
    s.reset()
    # store original
    if models[i].name() != 'STA':
        datalog = s.run(pulsetrain.characteristic_time(),
                        log=['isus.g','engine.time','membrane.V'],
                        log_interval=1)
        datalog = datalog.trim((nsteps)*period-tstep-20, (nsteps)*period, adjust=True)
        df = {'time': datalog['engine.time'],
              'gate': datalog.npview()['isus.g']/max(datalog['isus.g']),
              'sample': 0,
              'model': models[i].name(),
              'type': 'original'}
        df = pd.DataFrame(df)
        all_samples = all_samples.append(df, ignore_index=True)
    
    # re-calibrated
    for j, par_sample in enumerate(par_samples[i]):
        s.reset()
        for p, v in par_sample.items():
            name = p
            value = v
            if p.startswith("log"):
                name = p[4:]
                value = 10**v
            s.set_constant(name, value)
            
        # Log run time for comparisons
        t0 = time.time()
        datalog = s.run(pulsetrain.characteristic_time(),
                        log=['isus.g','engine.time','membrane.V'],
                        log_interval=1)
        t1 = time.time()
        dt = t1-t0
        
        datalog = datalog.trim((nsteps)*period-tstep-20, (nsteps)*period, adjust=True)
        df = {'time': datalog['engine.time'],
              'gate': datalog.npview()['isus.g']/max(datalog['isus.g']),
              'sample': j,
              'model': models[i].name(),
              'type': 'recalibrated',
              'runtime': dt}
        df = pd.DataFrame(df)
        all_samples = all_samples.append(df, ignore_index=True)



In [87]:

    
recalibrated_df = all_samples[all_samples['type']=='recalibrated']



In [88]:

    
sns.set(style="ticks")
sns.set_context('talk')

# Initialize the figure
f, ax = plt.subplots(figsize=(6, 3))

sns.boxplot(x='runtime', y='model', data=recalibrated_df,
            palette="pastel", whis="range")

# Tweak the visual presentation
ax.xaxis.grid(True)
ax.set(ylabel="")
ax.set(xlabel="time (ms)")
sns.despine(left=True)



In [89]:

    
#f.savefig('figures/isus/runtime_compare.pdf')



In [90]:

    
all_samples_detail = all_samples[(all_samples['time']>=19) & (all_samples['time'] < 50)]



In [91]:

    
sns.set_context('talk')

mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plt.subplots(figsize=(8,4))

g = sns.lineplot(x='time', y='gate', hue='model', style='type',
                 dashes = [(1,1),''],
                 data=all_samples, ax=ax,
                 legend=False)

ax2 = plt.axes([.6, .4, .2, .4])
sns.lineplot(x='time',y='gate',hue='model',style='type',
             dashes=[(1,1),''],
             data=all_samples_detail, ax=ax2,
             legend=False)
ax2.set_xlabel('')
ax2.set_ylabel('')

ax.set_ylabel('normalised gate')
ax.set_title('')









    Out[91]:





Text(0.5, 1.0, '')



In [92]:

    
#f.savefig('figures/isus/trace_compare.pdf')



In [ ]: