ABC calibration of $I_\text{CaL}$ in Courtemanche model to original dataset.


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 in the publication:

  • Steady-state activation [Li1997]
  • Steady-state inactivation [Li1997]
  • Inactivation time constant [Sun1997]

In [4]:
from experiments.ical_li import (li_act,
                                 li_inact_1000)
from experiments.ical_sun import sun_inact_kin_slow

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

Plot steady-state and time constant functions of original model


In [6]:
from ionchannelABC.visualization import plot_variables

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

V = np.arange(-140, 50, 0.01)

cou_par_map = {'di': 'ical.d_inf',
            'fi': 'ical.f_inf',
            'dt': 'ical.tau_d',
            'ft': 'ical.tau_f'}

f, ax = plot_variables(V, cou_par_map, 'models/courtemanche_ical.mmt', figshape=(2,2))


Activation gate ($d$) 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,
                                                li_act)

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

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


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 [12]:
limits = {'ical.p1': (-100, 100),
          'ical.p2': (0, 50),
          'log_ical.p3': (-7, 3),
          'ical.p4': (-100, 100),
          'ical.p5': (0, 50)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

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

Run ABC calibration


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

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

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


Theoretical minimum population size is 32 particles

In [17]:
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(2000),
             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.2446079782705148, '1': 1.2446079782705148, '2': 1.2446079782705148, '3': 1.2446079782705148, '4': 1.2446079782705148, '5': 1.2446079782705148, '6': 1.2446079782705148, '7': 1.042671915141306, '8': 0.6364361040472902, '9': 0.5199531035717904, '10': 0.40923240093228713, '11': 0.42429073603152656, '12': 0.521335957570653, '13': 1.2446079782705148, '14': 1.2446079782705148, '15': 1.2446079782705148}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=100, quantile_multiplier=1

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

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


INFO:History:Start <ABCSMC(id=1, start_time=2019-10-17 13:46:20.011046, end_time=None)>

In [20]:
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 2026
DEBUG:Epsilon:new eps, t=1, eps=1.5249660839874908
INFO:ABC:t:1 eps:1.5249660839874908
DEBUG:ABC:now submitting population 1
DEBUG:ABC:population 1 done
DEBUG:ABC:
total nr simulations up to t =1 is 5820
DEBUG:Epsilon:new eps, t=2, eps=0.8672636047163372
INFO:ABC:t:2 eps:0.8672636047163372
DEBUG:ABC:now submitting population 2
DEBUG:ABC:population 2 done
DEBUG:ABC:
total nr simulations up to t =2 is 9851
DEBUG:Epsilon:new eps, t=3, eps=0.6122752469319943
INFO:ABC:t:3 eps:0.6122752469319943
DEBUG:ABC:now submitting population 3
DEBUG:ABC:population 3 done
DEBUG:ABC:
total nr simulations up to t =3 is 14317
DEBUG:Epsilon:new eps, t=4, eps=0.5053136055782594
INFO:ABC:t:4 eps:0.5053136055782594
DEBUG:ABC:now submitting population 4
DEBUG:ABC:population 4 done
DEBUG:ABC:
total nr simulations up to t =4 is 19089
DEBUG:Epsilon:new eps, t=5, eps=0.43818344719701574
INFO:ABC:t:5 eps:0.43818344719701574
DEBUG:ABC:now submitting population 5
DEBUG:ABC:population 5 done
DEBUG:ABC:
total nr simulations up to t =5 is 24358
DEBUG:Epsilon:new eps, t=6, eps=0.35068839326566426
INFO:ABC:t:6 eps:0.35068839326566426
DEBUG:ABC:now submitting population 6
DEBUG:ABC:population 6 done
DEBUG:ABC:
total nr simulations up to t =6 is 29122
DEBUG:Epsilon:new eps, t=7, eps=0.26475756429715136
INFO:ABC:t:7 eps:0.26475756429715136
DEBUG:ABC:now submitting population 7
DEBUG:ABC:population 7 done
DEBUG:ABC:
total nr simulations up to t =7 is 33722
DEBUG:Epsilon:new eps, t=8, eps=0.1990459578923824
INFO:ABC:t:8 eps:0.1990459578923824
DEBUG:ABC:now submitting population 8
DEBUG:ABC:population 8 done
DEBUG:ABC:
total nr simulations up to t =8 is 38358
DEBUG:Epsilon:new eps, t=9, eps=0.14836609739749615
INFO:ABC:t:9 eps:0.14836609739749615
DEBUG:ABC:now submitting population 9
DEBUG:ABC:population 9 done
DEBUG:ABC:
total nr simulations up to t =9 is 43585
DEBUG:Epsilon:new eps, t=10, eps=0.1123759536913003
INFO:ABC:t:10 eps:0.1123759536913003
DEBUG:ABC:now submitting population 10
DEBUG:ABC:population 10 done
DEBUG:ABC:
total nr simulations up to t =10 is 48653
DEBUG:Epsilon:new eps, t=11, eps=0.08637742103956285
INFO:ABC:t:11 eps:0.08637742103956285
DEBUG:ABC:now submitting population 11
DEBUG:ABC:population 11 done
DEBUG:ABC:
total nr simulations up to t =11 is 53219
DEBUG:Epsilon:new eps, t=12, eps=0.06684543453822892
INFO:ABC:t:12 eps:0.06684543453822892
DEBUG:ABC:now submitting population 12
DEBUG:ABC:population 12 done
DEBUG:ABC:
total nr simulations up to t =12 is 64483
DEBUG:Epsilon:new eps, t=13, eps=0.05823023514550964
INFO:ABC:t:13 eps:0.05823023514550964
DEBUG:ABC:now submitting population 13
DEBUG:ABC:population 13 done
DEBUG:ABC:
total nr simulations up to t =13 is 70140
DEBUG:Epsilon:new eps, t=14, eps=0.05363456607282805
INFO:ABC:t:14 eps:0.05363456607282805
DEBUG:ABC:now submitting population 14
DEBUG:ABC:population 14 done
DEBUG:ABC:
total nr simulations up to t =14 is 83394
DEBUG:Epsilon:new eps, t=15, eps=0.050600813809692904
INFO:ABC:t:15 eps:0.050600813809692904
DEBUG:ABC:now submitting population 15
DEBUG:ABC:population 15 done
DEBUG:ABC:
total nr simulations up to t =15 is 91162
DEBUG:Epsilon:new eps, t=16, eps=0.04893068140780184
INFO:ABC:t:16 eps:0.04893068140780184
DEBUG:ABC:now submitting population 16
DEBUG:ABC:population 16 done
DEBUG:ABC:
total nr simulations up to t =16 is 105966
DEBUG:Epsilon:new eps, t=17, eps=0.04801433437681319
INFO:ABC:t:17 eps:0.04801433437681319
DEBUG:ABC:now submitting population 17
DEBUG:ABC:population 17 done
DEBUG:ABC:
total nr simulations up to t =17 is 133124
DEBUG:Epsilon:new eps, t=18, eps=0.04737246637920021
INFO:ABC:t:18 eps:0.04737246637920021
DEBUG:ABC:now submitting population 18
DEBUG:ABC:population 18 done
DEBUG:ABC:
total nr simulations up to t =18 is 254528
DEBUG:Epsilon:new eps, t=19, eps=0.043859848176284015
INFO:ABC:t:19 eps:0.043859848176284015
DEBUG:ABC:now submitting population 19
DEBUG:ABC:population 19 done
DEBUG:ABC:
total nr simulations up to t =19 is 434651
DEBUG:Epsilon:new eps, t=20, eps=0.03911542472357917
INFO:ABC:t:20 eps:0.03911542472357917
DEBUG:ABC:now submitting population 20
DEBUG:ABC:population 20 done
DEBUG:ABC:
total nr simulations up to t =20 is 509945
DEBUG:Epsilon:new eps, t=21, eps=0.03563553543525459
INFO:ABC:t:21 eps:0.03563553543525459
DEBUG:ABC:now submitting population 21
DEBUG:ABC:population 21 done
DEBUG:ABC:
total nr simulations up to t =21 is 600128
DEBUG:Epsilon:new eps, t=22, eps=0.0328135014894609
INFO:ABC:t:22 eps:0.0328135014894609
DEBUG:ABC:now submitting population 22
DEBUG:ABC:population 22 done
DEBUG:ABC:
total nr simulations up to t =22 is 726072
DEBUG:Epsilon:new eps, t=23, eps=0.030547365048097404
INFO:ABC:t:23 eps:0.030547365048097404
DEBUG:ABC:now submitting population 23
DEBUG:ABC:population 23 done
DEBUG:ABC:
total nr simulations up to t =23 is 895877
DEBUG:Epsilon:new eps, t=24, eps=0.02892550724789828
INFO:ABC:t:24 eps:0.02892550724789828
DEBUG:ABC:now submitting population 24
DEBUG:ABC:population 24 done
DEBUG:ABC:
total nr simulations up to t =24 is 1079087
DEBUG:Epsilon:new eps, t=25, eps=0.02779584784530078
INFO:ABC:t:25 eps:0.02779584784530078
DEBUG:ABC:now submitting population 25
DEBUG:ABC:population 25 done
DEBUG:ABC:
total nr simulations up to t =25 is 1420182
DEBUG:Epsilon:new eps, t=26, eps=0.026807807548790873
INFO:History:Done <ABCSMC(id=1, start_time=2019-10-17 13:46:20.011046, end_time=2019-10-17 15:10:22.909217)>

Analysis of results


In [17]:
history = History('sqlite:///results/courtemanche/ical/original/courtemanche_ical_dgate_original.db')

In [18]:
df, w = history.get_distribution(m=0)

In [19]:
df.describe()


Out[19]:
name ical.p1 ical.p2 ical.p4 ical.p5 log_ical.p3
count 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000
mean -64.968145 10.380238 17.914814 6.142103 -6.351218
std 15.102658 5.168971 4.096174 0.427570 0.969497
min -99.203878 0.228430 3.928299 5.615548 -6.807049
25% -75.280129 6.162337 18.678294 5.916643 -6.730432
50% -66.906606 10.211761 19.252756 6.037111 -6.696019
75% -56.688353 14.655124 19.780059 6.182823 -6.648474
max -19.419167 22.396021 21.104552 8.742910 -2.035612

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

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

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

plt.tight_layout()



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

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

f, ax = plot_variables(V, cou_par_map, 
                       'models/courtemanche_ical.mmt', 
                       [cou_par_samples],
                       figshape=(2,2))
plt.tight_layout()



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

In [24]:
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 [25]:
sns.set_context('paper')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)
plt.tight_layout()


Voltage-dependent inactivation gate ($f$) calibration


In [26]:
observations, model, summary_statistics = setup(modelfile,
                                                li_inact_1000,
                                                sun_inact_kin_slow)

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

In [28]:
g = plot_sim_results(modelfile,
                     li_inact_1000,
                     sun_inact_kin_slow)


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 [29]:
limits = {'log_ical.q1': (0, 3),
          'log_ical.q2': (-2, 3),
          'log_ical.q3': (-4, 0),
          'ical.q4': (-100, 100),
          'log_ical.q5': (-4, 0),
          'ical.q6': (-100, 100),
          'ical.q7': (0, 50)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

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

Run ABC calibration


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

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

In [36]:
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 [37]:
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(2000),
             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.176132058445145, '1': 1.176132058445145, '2': 1.176132058445145, '3': 1.176132058445145, '4': 0.3892784268190627, '5': 0.604405978482229, '6': 0.8074486118786046, '7': 1.0034312846646722, '8': 1.176132058445145, '9': 1.176132058445145, '10': 1.176132058445145, '11': 1.0995044927708642, '12': 1.0034312846646722, '13': 0.7544045424851181, '14': 0.5263305611568131, '15': 1.104247498336103, '16': 1.1063428446327512, '17': 1.8421569640488558, '18': 0.5260931009442432}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=100, quantile_multiplier=1

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

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


INFO:History:Start <ABCSMC(id=1, start_time=2019-10-17 15:25:47.936468, 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 10309
DEBUG:Epsilon:new eps, t=1, eps=8.705169594851025
INFO:ABC:t:1 eps:8.705169594851025
DEBUG:ABC:now submitting population 1
DEBUG:ABC:population 1 done
DEBUG:ABC:
total nr simulations up to t =1 is 15452
DEBUG:Epsilon:new eps, t=2, eps=2.1686344483401947
INFO:ABC:t:2 eps:2.1686344483401947
DEBUG:ABC:now submitting population 2
DEBUG:ABC:population 2 done
DEBUG:ABC:
total nr simulations up to t =2 is 20832
DEBUG:Epsilon:new eps, t=3, eps=1.4965119735542094
INFO:ABC:t:3 eps:1.4965119735542094
DEBUG:ABC:now submitting population 3
DEBUG:ABC:population 3 done
DEBUG:ABC:
total nr simulations up to t =3 is 28119
DEBUG:Epsilon:new eps, t=4, eps=1.153585553861291
INFO:ABC:t:4 eps:1.153585553861291
DEBUG:ABC:now submitting population 4
DEBUG:ABC:population 4 done
DEBUG:ABC:
total nr simulations up to t =4 is 36782
DEBUG:Epsilon:new eps, t=5, eps=0.8765080578815893
INFO:ABC:t:5 eps:0.8765080578815893
DEBUG:ABC:now submitting population 5
DEBUG:ABC:population 5 done
DEBUG:ABC:
total nr simulations up to t =5 is 48795
DEBUG:Epsilon:new eps, t=6, eps=0.691020560153727
INFO:ABC:t:6 eps:0.691020560153727
DEBUG:ABC:now submitting population 6
DEBUG:ABC:population 6 done
DEBUG:ABC:
total nr simulations up to t =6 is 64716
DEBUG:Epsilon:new eps, t=7, eps=0.5650043476709863
INFO:ABC:t:7 eps:0.5650043476709863
DEBUG:ABC:now submitting population 7
DEBUG:ABC:population 7 done
DEBUG:ABC:
total nr simulations up to t =7 is 83632
DEBUG:Epsilon:new eps, t=8, eps=0.47680147302527776
INFO:ABC:t:8 eps:0.47680147302527776
DEBUG:ABC:now submitting population 8

Analysis of results


In [31]:
history = History('sqlite:///results/courtemanche/ical/original/courtemanche_ical_fgate_original.db')

In [32]:
history.all_runs()


Out[32]:
[<ABCSMC(id=1, start_time=2019-10-17 15:25:47.936468, end_time=2019-10-17 19:38:25.461218)>]

In [33]:
df, w = history.get_distribution(m=0)

In [34]:
df.describe()


Out[34]:
name ical.q4 ical.q6 ical.q7 log_ical.q1 log_ical.q2 log_ical.q3 log_ical.q5
count 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000
mean -15.703341 27.702654 7.340350 1.590756 -0.481084 -2.109852 -1.782001
std 35.742192 1.485608 1.000522 0.710907 0.847052 0.803497 1.038426
min -99.905497 22.728886 3.534692 0.028940 -1.989594 -3.995052 -3.996133
25% -32.662145 26.862679 6.779595 1.047811 -1.141645 -2.601428 -2.579144
50% -14.199977 27.544204 7.320581 1.630087 -0.499407 -1.947088 -1.666768
75% -5.441050 28.246457 7.836125 2.125888 0.123311 -1.600957 -0.935848
max 99.841123 35.207710 12.626933 2.996731 2.871118 -0.007737 -0.004483

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

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

g = plot_sim_results(modelfile,
                     li_inact_1000,
                     sun_inact_kin_slow,
                     df=df, w=w)

plt.tight_layout()



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 [39]:
sns.set_context('paper')
g = plot_kde_matrix_custom(df, w, limits=limits, refval=originals)
plt.tight_layout()



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

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

f, ax = plot_variables(V, cou_par_map, 
                       'models/courtemanche_ical.mmt', 
                       [cou_par_samples],
                       figshape=(2,2))
plt.tight_layout()