ABC calibration of $I_\text{Kur}$ in Nygren model to original dataset.

Note the term $I_\text{sus}$ for sustained outward Potassium current is used throughout the notebook.


In [1]:
import os, tempfile
import logging
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

In [2]:
from ionchannelABC import theoretical_population_size
from ionchannelABC import IonChannelDistance, EfficientMultivariateNormalTransition, IonChannelAcceptor
from ionchannelABC.experiment import setup
from ionchannelABC.visualization import plot_sim_results, plot_kde_matrix_custom
import myokit


INFO:myokit:Loading Myokit version 1.28.3

In [3]:
from pyabc import Distribution, RV, History, ABCSMC
from pyabc.epsilon import MedianEpsilon
from pyabc.sampler import MulticoreEvalParallelSampler, SingleCoreSampler
from pyabc.populationstrategy import ConstantPopulationSize

Initial set-up

Load experiments used for original dataset calibration:

  • Steady-state activation [Wang1993]
  • Activation time constant [Wang1993]
  • Steady-state inactivation [Firek1995]
  • Inactivation time constant [Nygren1998]
  • Recovery time constant [Nygren1998]

In [4]:
from experiments.isus_wang import wang_act_and_kin
from experiments.isus_firek import (firek_inact)
from experiments.isus_nygren import (nygren_inact_kin,
                                     nygren_rec)

In [5]:
modelfile = 'models/nygren_isus.mmt'

Plot steady-state and time constant functions of original model


In [6]:
from ionchannelABC.visualization import plot_variables

In [7]:
sns.set_context('talk')

V = np.arange(-100, 40, 0.01)

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

f, ax = plot_variables(V, nyg_par_map, modelfile, figshape=(2,2))


Activation gate ($r$) calibration

Combine model and experiments to produce:

  • observations dataframe
  • model function to run experiments and return traces
  • summary statistics function to accept traces

In [8]:
observations, model, summary_statistics = setup(modelfile,
                                                wang_act_and_kin)

In [9]:
assert len(observations)==len(summary_statistics(model({})))

In [10]:
g = plot_sim_results(modelfile,
                     wang_act_and_kin)


Set up prior ranges for each parameter in the model.

See the modelfile for further information on specific parameters. Prepending `log_' has the effect of setting the parameter in log space.


In [11]:
limits = {'isus.p1': (-100, 100),
          'isus.p2': (1e-7, 50),
          'log_isus.p3': (-5, 0),
          'isus.p4': (-100, 100),
          'isus.p5': (1e-7, 50),
          'log_isus.p6': (-6, -1)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

In [12]:
# Test this works correctly with set-up functions
assert len(observations) == len(summary_statistics(model(prior.rvs())))

Run ABC calibration


In [12]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "nygren_isus_rgate_original.db"))

In [13]:
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)

In [14]:
pop_size = theoretical_population_size(2, len(limits))
print("Theoretical minimum population size is {} particles".format(pop_size))


Theoretical minimum population size is 64 particles

In [15]:
abc = ABCSMC(models=model,
             parameter_priors=prior,
             distance_function=IonChannelDistance(
                 exp_id=list(observations.exp_id),
                 variance=list(observations.variance),
                 delta=0.05),
             population_size=ConstantPopulationSize(500),
             summary_statistics=summary_statistics,
             transitions=EfficientMultivariateNormalTransition(),
             eps=MedianEpsilon(initial_epsilon=100),
             sampler=MulticoreEvalParallelSampler(n_procs=16),
             acceptor=IonChannelAcceptor())


DEBUG:ABC:ion channel weights: {'0': 1.2435658718644211, '1': 1.243311823678135, '2': 1.5188448657111304, '3': 1.519224008213791, '4': 1.3670286111645151, '5': 0.8044958859093764, '6': 1.1397025050382832, '7': 1.1399159724488552, '8': 1.1399159724488594, '9': 1.1397025050382832, '10': 0.08446935501503675, '11': 0.32315877006716853, '12': 0.3912264863854337, '13': 0.8745062636850891, '14': 1.1442949744557158, '15': 0.9908407438044368, '16': 0.8745062636850831, '17': 1.061289121386388}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=100, quantile_multiplier=1

In [16]:
obs = observations.to_dict()['y']
obs = {str(k): v for k, v in obs.items()}

In [17]:
abc_id = abc.new(db_path, obs)


INFO:History:Start <ABCSMC(id=1, start_time=2019-10-30 09:52:42.170175, end_time=None)>

In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)


INFO:ABC:t:0 eps:100
DEBUG:ABC:now submitting population 0
DEBUG:ABC:population 0 done
DEBUG:ABC:
total nr simulations up to t =0 is 948
DEBUG:Epsilon:new eps, t=1, eps=1.8482362538424795
INFO:ABC:t:1 eps:1.8482362538424795
DEBUG:ABC:now submitting population 1
DEBUG:ABC:population 1 done
DEBUG:ABC:
total nr simulations up to t =1 is 2022
DEBUG:Epsilon:new eps, t=2, eps=1.2559647857588863
INFO:ABC:t:2 eps:1.2559647857588863
DEBUG:ABC:now submitting population 2
DEBUG:ABC:population 2 done
DEBUG:ABC:
total nr simulations up to t =2 is 3185
DEBUG:Epsilon:new eps, t=3, eps=1.0877007372731602
INFO:ABC:t:3 eps:1.0877007372731602
DEBUG:ABC:now submitting population 3
DEBUG:ABC:population 3 done
DEBUG:ABC:
total nr simulations up to t =3 is 4422
DEBUG:Epsilon:new eps, t=4, eps=1.0338740966404782
INFO:ABC:t:4 eps:1.0338740966404782
DEBUG:ABC:now submitting population 4

Analysis of results


In [13]:
history = History('sqlite:///results/nygren/isus/original/nygren_isus_rgate_original.db')

In [14]:
history.all_runs()


Out[14]:
[<ABCSMC(id=1, start_time=2019-10-30 09:52:42.170175, end_time=2019-10-30 10:51:45.325455)>]

In [15]:
df, w = history.get_distribution()

In [16]:
df.describe()


Out[16]:
name isus.p1 isus.p2 isus.p4 isus.p5 log_isus.p3 log_isus.p6
count 500.000000 500.000000 500.000000 500.000000 500.000000 500.000000
mean 1.891587 4.543242 -13.606205 7.628083 -2.384896 -3.107227
std 0.031980 0.020526 0.313763 0.200001 0.004245 0.018890
min 1.816604 4.415867 -14.425018 7.059888 -2.396415 -3.155123
25% 1.869075 4.536539 -13.838517 7.487667 -2.387434 -3.120984
50% 1.890769 4.547424 -13.614881 7.638881 -2.384380 -3.108035
75% 1.912584 4.556689 -13.409248 7.772041 -2.381918 -3.094375
max 1.979269 4.575717 -12.821089 8.179115 -2.374562 -3.060003

In [17]:
sns.set_context('poster')

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

g = plot_sim_results(modelfile,
                     wang_act_and_kin,
                     df=df, w=w)

plt.tight_layout()



In [18]:
import pandas as pd
N = 100
nyg_par_samples = df.sample(n=N, weights=w, replace=True)
nyg_par_samples = nyg_par_samples.set_index([pd.Index(range(N))])
nyg_par_samples = nyg_par_samples.to_dict(orient='records')

In [19]:
sns.set_context('talk')
mpl.rcParams['font.size'] = 14
mpl.rcParams['legend.fontsize'] = 14

f, ax = plot_variables(V, nyg_par_map, 
                       'models/nygren_isus.mmt', 
                       [nyg_par_samples],
                       figshape=(2,2))



In [20]:
from ionchannelABC.visualization import plot_kde_matrix_custom
import myokit
import numpy as np

In [21]:
m,_,_ = myokit.load(modelfile)

In [22]:
originals = {}
for name in limits.keys():
    if name.startswith("log"):
        name_ = name[4:]
    else:
        name_ = name
    val = m.value(name_)
    if name.startswith("log"):
        val_ = np.log10(val)
    else:
        val_ = val
    originals[name] = val_

In [23]:
sns.set_context('paper')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)
plt.tight_layout()


Inactivation gate ($s$) calibration


In [24]:
observations, model, summary_statistics = setup(modelfile,
                                                firek_inact,
                                                nygren_inact_kin,
                                                nygren_rec)

In [25]:
assert len(observations)==len(summary_statistics(model({})))

In [27]:
sns.set_context('talk')
g = plot_sim_results(modelfile,
                     firek_inact,
                     nygren_inact_kin,
                     nygren_rec)



In [28]:
limits = {'isus.q1': (-100, 100),
          'isus.q2': (1e-7, 50),
          'isus.q3': (0., 1.),
          'log_isus.q4': (-4, 1),
          'isus.q5': (-100, 100),
          'isus.q6': (1e-7, 50),
          'log_isus.q7': (-3, 2)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

In [29]:
# Test this works correctly with set-up functions
assert len(observations) == len(summary_statistics(model(prior.rvs())))

Run ABC calibration


In [22]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "nygren_isus_sgate_original.db"))

In [23]:
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)

In [24]:
pop_size = theoretical_population_size(2, len(limits))
print("Theoretical minimum population size is {} particles".format(pop_size))


Theoretical minimum population size is 128 particles

In [18]:
abc = ABCSMC(models=model,
             parameter_priors=prior,
             distance_function=IonChannelDistance(
                 exp_id=list(observations.exp_id),
                 variance=list(observations.variance),
                 delta=0.05),
             population_size=ConstantPopulationSize(1000),
             summary_statistics=summary_statistics,
             transitions=EfficientMultivariateNormalTransition(),
             eps=MedianEpsilon(initial_epsilon=100),
             sampler=MulticoreEvalParallelSampler(n_procs=16),
             acceptor=IonChannelAcceptor())


DEBUG:ABC:ion channel weights: {'0': 0.4248110530506898, '1': 0.9476554260361589, '2': 0.47382771301807824, '3': 0.513313355769584, '4': 0.35198630109914286, '5': 0.6159760269235005, '6': 0.947655426036154, '7': 1.0466865780553765, '8': 1.1324846983182664, '9': 0.6459809613606208, '10': 1.0541409114085802, '11': 0.9789135226692829, '12': 0.9906822580908385, '13': 3.875885768163725}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=100, quantile_multiplier=1

In [19]:
obs = observations.to_dict()['y']
obs = {str(k): v for k, v in obs.items()}

In [20]:
abc_id = abc.new(db_path, obs)


INFO:History:Start <ABCSMC(id=4, start_time=2019-10-28 15:25:42.328733, end_time=None)>

In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)


INFO:ABC:t:0 eps:100
DEBUG:ABC:now submitting population 0

Database results analysis


In [30]:
history = History('sqlite:///results/nygren/isus/original/nygren_isus_sgate_original.db')

In [31]:
df, w = history.get_distribution()

In [32]:
df.describe()


Out[32]:
name isus.q1 isus.q2 isus.q3 isus.q5 isus.q6 log_isus.q4 log_isus.q7
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 24.236605 13.347338 0.139650 41.453206 24.475241 -1.176228 -0.568873
std 11.967437 7.079625 0.083050 44.280695 13.780677 0.238996 0.121901
min 5.080157 0.161439 0.000073 -99.113497 0.032287 -1.837425 -2.512781
25% 16.538815 8.110749 0.069655 20.057140 12.294586 -1.334780 -0.565752
50% 21.464245 12.504027 0.134244 49.420872 24.216142 -1.186948 -0.542518
75% 27.973434 17.571098 0.201953 74.191748 36.410594 -1.032830 -0.525306
max 89.267139 36.818964 0.342169 99.987830 49.988193 -0.398952 -0.495444

In [33]:
sns.set_context('poster')

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

g = plot_sim_results(modelfile,
                     firek_inact,
                     nygren_inact_kin,
                     nygren_rec,
                     df=df, w=w)

plt.tight_layout()



In [34]:
N = 100
nyg_par_samples = df.sample(n=N, weights=w, replace=True)
nyg_par_samples = nyg_par_samples.set_index([pd.Index(range(N))])
nyg_par_samples = nyg_par_samples.to_dict(orient='records')

In [35]:
sns.set_context('talk')

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

f, ax = plot_variables(V, nyg_par_map, 
                       'models/nygren_isus.mmt', 
                       [nyg_par_samples],
                       figshape=(2,2))



In [36]:
m,_,_ = myokit.load(modelfile)

In [37]:
originals = {}
for name in limits.keys():
    if name.startswith("log"):
        name_ = name[4:]
    else:
        name_ = name
    val = m.value(name_)
    if name.startswith("log"):
        val_ = np.log10(val)
    else:
        val_ = val
    originals[name] = val_

In [38]:
sns.set_context('notebook')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)



In [ ]: