ABC calibration of $I_\text{Na}$ 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.29.1

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 [Sakakibara1992]
  • Activation time constant [Schneider1994]
  • Steady-state inactivation [Sakakibara1992]
  • Inactivation time constant [Sakakibara1992]
  • Recovery time constant [Sakakibara1992]

In [4]:
from experiments.ina_sakakibara import (sakakibara_act,
                                        sakakibara_inact,
                                        sakakibara_inact_kin_fast,
                                        sakakibara_inact_kin_slow,
                                        sakakibara_rec_fast,
                                        sakakibara_rec_slow)
from experiments.ina_schneider import (schneider_taum)

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

Activation gate ($m$) calibration

Combine model and experiments to produce:

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

In [16]:
observations, model, summary_statistics = setup(modelfile,
                                                sakakibara_act,
                                                schneider_taum_cou_adjust)

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

In [19]:
# Test the output of the unaltered model.
g = plot_sim_results(modelfile,
                     sakakibara_act,
                     schneider_taum_cou_adjust)


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 [24]:
limits = {'ina.a1_m': (-100, 0),
          'ina.a2_m': (0, 1),
          'ina.a3_m': (0, 1),
          'ina.a4_m': (0, 10),
          'ina.b1_m': (0, 10),
          'ina.b2_m': (0, 100)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

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

Run ABC calibration


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

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

Test theoretical number of particles for approximately 2 particles per dimension in the initial sampling of the parameter hyperspace.


In [28]:
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

Initialise ABCSMC (see pyABC documentation for further details).

IonChannelDistance calculates the weighting applied to each datapoint based on the experimental variance.


In [30]:
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=20),
             sampler=MulticoreEvalParallelSampler(n_procs=16),
             acceptor=IonChannelAcceptor())


DEBUG:ABC:ion channel weights: {'0': 1.5570757359362264, '1': 1.5570757359362264, '2': 1.5570757359362264, '3': 1.5570757359362264, '4': 1.5570757359362264, '5': 0.5904911622959486, '6': 0.45939980914873735, '7': 0.5710013425594168, '8': 1.5570757359362264, '9': 1.5570757359362264, '10': 1.5570757359362264, '11': 1.5570757359362264, '12': 1.5570757359362264, '13': 0.15235445469736852, '14': 0.1819349064747563, '15': 0.1879980741832115, '16': 0.25507625823484376, '17': 0.4491978085688121, '18': 0.5936989296806218, '19': 1.3185718986681156, '20': 0.4160669679218102, '21': 1.2534510282040947}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=20, quantile_multiplier=1

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

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


INFO:History:Start <ABCSMC(id=1, start_time=2019-12-19 11:31:55.095938, end_time=None)>

Run calibration with stopping criterion of particle 1\% acceptance rate.


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


INFO:ABC:t: 0, eps: 20.
DEBUG:ABC:Now submitting population 0.
DEBUG:ABC:Population 0 done.
DEBUG:ABC:Total samples up to t = 0: 3517.
INFO:ABC:Acceptance rate: 2000 / 3517 = 5.6867e-01.
DEBUG:Epsilon:new eps, t=1, eps=1.1920110107351267
INFO:ABC:t: 1, eps: 1.1920110107351267.
DEBUG:ABC:Now submitting population 1.
DEBUG:ABC:Population 1 done.
DEBUG:ABC:Total samples up to t = 1: 7776.
INFO:ABC:Acceptance rate: 2000 / 4259 = 4.6959e-01.
DEBUG:Epsilon:new eps, t=2, eps=0.8990569694575163
INFO:ABC:t: 2, eps: 0.8990569694575163.
DEBUG:ABC:Now submitting population 2.
DEBUG:ABC:Population 2 done.
DEBUG:ABC:Total samples up to t = 2: 12190.
INFO:ABC:Acceptance rate: 2000 / 4414 = 4.5310e-01.
DEBUG:Epsilon:new eps, t=3, eps=0.8096494247442982
INFO:ABC:t: 3, eps: 0.8096494247442982.
DEBUG:ABC:Now submitting population 3.
DEBUG:ABC:Population 3 done.
DEBUG:ABC:Total samples up to t = 3: 17388.
INFO:ABC:Acceptance rate: 2000 / 5198 = 3.8476e-01.
DEBUG:Epsilon:new eps, t=4, eps=0.7508559436975845
INFO:ABC:t: 4, eps: 0.7508559436975845.
DEBUG:ABC:Now submitting population 4.
DEBUG:ABC:Population 4 done.
DEBUG:ABC:Total samples up to t = 4: 22684.
INFO:ABC:Acceptance rate: 2000 / 5296 = 3.7764e-01.
DEBUG:Epsilon:new eps, t=5, eps=0.7041252958957414
INFO:ABC:t: 5, eps: 0.7041252958957414.
DEBUG:ABC:Now submitting population 5.
DEBUG:ABC:Population 5 done.
DEBUG:ABC:Total samples up to t = 5: 28590.
INFO:ABC:Acceptance rate: 2000 / 5906 = 3.3864e-01.
DEBUG:Epsilon:new eps, t=6, eps=0.662069237730507
INFO:ABC:t: 6, eps: 0.662069237730507.
DEBUG:ABC:Now submitting population 6.
DEBUG:ABC:Population 6 done.
DEBUG:ABC:Total samples up to t = 6: 34571.
INFO:ABC:Acceptance rate: 2000 / 5981 = 3.3439e-01.
DEBUG:Epsilon:new eps, t=7, eps=0.622661070270785
INFO:ABC:t: 7, eps: 0.622661070270785.
DEBUG:ABC:Now submitting population 7.
DEBUG:ABC:Population 7 done.
DEBUG:ABC:Total samples up to t = 7: 41055.
INFO:ABC:Acceptance rate: 2000 / 6484 = 3.0845e-01.
DEBUG:Epsilon:new eps, t=8, eps=0.5909612165818797
INFO:ABC:t: 8, eps: 0.5909612165818797.
DEBUG:ABC:Now submitting population 8.
DEBUG:ABC:Population 8 done.
DEBUG:ABC:Total samples up to t = 8: 47674.
INFO:ABC:Acceptance rate: 2000 / 6619 = 3.0216e-01.
DEBUG:Epsilon:new eps, t=9, eps=0.5698579088020338
INFO:ABC:t: 9, eps: 0.5698579088020338.
DEBUG:ABC:Now submitting population 9.
DEBUG:ABC:Population 9 done.
DEBUG:ABC:Total samples up to t = 9: 56841.
INFO:ABC:Acceptance rate: 2000 / 9167 = 2.1817e-01.
DEBUG:Epsilon:new eps, t=10, eps=0.5481013686200193
INFO:ABC:t: 10, eps: 0.5481013686200193.
DEBUG:ABC:Now submitting population 10.
DEBUG:ABC:Population 10 done.
DEBUG:ABC:Total samples up to t = 10: 64524.
INFO:ABC:Acceptance rate: 2000 / 7683 = 2.6031e-01.
DEBUG:Epsilon:new eps, t=11, eps=0.5310636150592549
INFO:ABC:t: 11, eps: 0.5310636150592549.
DEBUG:ABC:Now submitting population 11.
DEBUG:ABC:Population 11 done.
DEBUG:ABC:Total samples up to t = 11: 72756.
INFO:ABC:Acceptance rate: 2000 / 8232 = 2.4295e-01.
DEBUG:Epsilon:new eps, t=12, eps=0.5155487795509347
INFO:ABC:t: 12, eps: 0.5155487795509347.
DEBUG:ABC:Now submitting population 12.
DEBUG:ABC:Population 12 done.
DEBUG:ABC:Total samples up to t = 12: 81954.
INFO:ABC:Acceptance rate: 2000 / 9198 = 2.1744e-01.
DEBUG:Epsilon:new eps, t=13, eps=0.5028328179256911
INFO:ABC:t: 13, eps: 0.5028328179256911.
DEBUG:ABC:Now submitting population 13.
DEBUG:ABC:Population 13 done.
DEBUG:ABC:Total samples up to t = 13: 90430.
INFO:ABC:Acceptance rate: 2000 / 8476 = 2.3596e-01.
DEBUG:Epsilon:new eps, t=14, eps=0.4923246418918639
INFO:ABC:t: 14, eps: 0.4923246418918639.
DEBUG:ABC:Now submitting population 14.
DEBUG:ABC:Population 14 done.
DEBUG:ABC:Total samples up to t = 14: 99014.
INFO:ABC:Acceptance rate: 2000 / 8584 = 2.3299e-01.
DEBUG:Epsilon:new eps, t=15, eps=0.4830946464943411
INFO:ABC:t: 15, eps: 0.4830946464943411.
DEBUG:ABC:Now submitting population 15.
DEBUG:ABC:Population 15 done.
DEBUG:ABC:Total samples up to t = 15: 107899.
INFO:ABC:Acceptance rate: 2000 / 8885 = 2.2510e-01.
DEBUG:Epsilon:new eps, t=16, eps=0.47540657387160645
INFO:ABC:t: 16, eps: 0.47540657387160645.
DEBUG:ABC:Now submitting population 16.
DEBUG:ABC:Population 16 done.
DEBUG:ABC:Total samples up to t = 16: 119297.
INFO:ABC:Acceptance rate: 2000 / 11398 = 1.7547e-01.
DEBUG:Epsilon:new eps, t=17, eps=0.4687071617589745
INFO:ABC:t: 17, eps: 0.4687071617589745.
DEBUG:ABC:Now submitting population 17.
DEBUG:ABC:Population 17 done.
DEBUG:ABC:Total samples up to t = 17: 133750.
INFO:ABC:Acceptance rate: 2000 / 14453 = 1.3838e-01.
DEBUG:Epsilon:new eps, t=18, eps=0.46333993353289404
INFO:ABC:t: 18, eps: 0.46333993353289404.
DEBUG:ABC:Now submitting population 18.
DEBUG:ABC:Population 18 done.
DEBUG:ABC:Total samples up to t = 18: 153126.
INFO:ABC:Acceptance rate: 2000 / 19376 = 1.0322e-01.
DEBUG:Epsilon:new eps, t=19, eps=0.45799442399995893
INFO:ABC:t: 19, eps: 0.45799442399995893.
DEBUG:ABC:Now submitting population 19.
DEBUG:ABC:Population 19 done.
DEBUG:ABC:Total samples up to t = 19: 182852.
INFO:ABC:Acceptance rate: 2000 / 29726 = 6.7281e-02.
DEBUG:Epsilon:new eps, t=20, eps=0.4529371174398475
INFO:ABC:t: 20, eps: 0.4529371174398475.
DEBUG:ABC:Now submitting population 20.
DEBUG:ABC:Population 20 done.
DEBUG:ABC:Total samples up to t = 20: 225670.
INFO:ABC:Acceptance rate: 2000 / 42818 = 4.6709e-02.
DEBUG:Epsilon:new eps, t=21, eps=0.4486546449731911
INFO:ABC:t: 21, eps: 0.4486546449731911.
DEBUG:ABC:Now submitting population 21.
DEBUG:ABC:Population 21 done.
DEBUG:ABC:Total samples up to t = 21: 272267.
INFO:ABC:Acceptance rate: 2000 / 46597 = 4.2921e-02.
DEBUG:Epsilon:new eps, t=22, eps=0.44522588242339795
INFO:ABC:t: 22, eps: 0.44522588242339795.
DEBUG:ABC:Now submitting population 22.
DEBUG:ABC:Population 22 done.
DEBUG:ABC:Total samples up to t = 22: 318255.
INFO:ABC:Acceptance rate: 2000 / 45988 = 4.3490e-02.
DEBUG:Epsilon:new eps, t=23, eps=0.44173058274877325
INFO:ABC:t: 23, eps: 0.44173058274877325.
DEBUG:ABC:Now submitting population 23.
DEBUG:ABC:Population 23 done.
DEBUG:ABC:Total samples up to t = 23: 373837.
INFO:ABC:Acceptance rate: 2000 / 55582 = 3.5983e-02.
DEBUG:Epsilon:new eps, t=24, eps=0.4387899435617417
INFO:ABC:t: 24, eps: 0.4387899435617417.
DEBUG:ABC:Now submitting population 24.
DEBUG:ABC:Population 24 done.
DEBUG:ABC:Total samples up to t = 24: 439976.
INFO:ABC:Acceptance rate: 2000 / 66139 = 3.0239e-02.
DEBUG:Epsilon:new eps, t=25, eps=0.4363575709273443
INFO:ABC:t: 25, eps: 0.4363575709273443.
DEBUG:ABC:Now submitting population 25.
DEBUG:ABC:Population 25 done.
DEBUG:ABC:Total samples up to t = 25: 516210.
INFO:ABC:Acceptance rate: 2000 / 76234 = 2.6235e-02.
DEBUG:Epsilon:new eps, t=26, eps=0.4342825704887904
INFO:ABC:t: 26, eps: 0.4342825704887904.
DEBUG:ABC:Now submitting population 26.
DEBUG:ABC:Population 26 done.
DEBUG:ABC:Total samples up to t = 26: 598385.
INFO:ABC:Acceptance rate: 2000 / 82175 = 2.4338e-02.
DEBUG:Epsilon:new eps, t=27, eps=0.4325727136384069
INFO:ABC:t: 27, eps: 0.4325727136384069.
DEBUG:ABC:Now submitting population 27.
DEBUG:ABC:Population 27 done.
DEBUG:ABC:Total samples up to t = 27: 673833.
INFO:ABC:Acceptance rate: 2000 / 75448 = 2.6508e-02.
DEBUG:Epsilon:new eps, t=28, eps=0.4311006269116266
INFO:ABC:t: 28, eps: 0.4311006269116266.
DEBUG:ABC:Now submitting population 28.
DEBUG:ABC:Population 28 done.
DEBUG:ABC:Total samples up to t = 28: 732772.
INFO:ABC:Acceptance rate: 2000 / 58939 = 3.3933e-02.
DEBUG:Epsilon:new eps, t=29, eps=0.4298413772185078
INFO:ABC:t: 29, eps: 0.4298413772185078.
DEBUG:ABC:Now submitting population 29.
DEBUG:ABC:Population 29 done.
DEBUG:ABC:Total samples up to t = 29: 774386.
INFO:ABC:Acceptance rate: 2000 / 41614 = 4.8061e-02.
DEBUG:Epsilon:new eps, t=30, eps=0.42888493600174343
INFO:ABC:t: 30, eps: 0.42888493600174343.
DEBUG:ABC:Now submitting population 30.
DEBUG:ABC:Population 30 done.
DEBUG:ABC:Total samples up to t = 30: 817308.
INFO:ABC:Acceptance rate: 2000 / 42922 = 4.6596e-02.
DEBUG:Epsilon:new eps, t=31, eps=0.4281682830995656
INFO:ABC:t: 31, eps: 0.4281682830995656.
DEBUG:ABC:Now submitting population 31.
DEBUG:ABC:Population 31 done.
DEBUG:ABC:Total samples up to t = 31: 828955.
INFO:ABC:Acceptance rate: 2000 / 11647 = 1.7172e-01.
DEBUG:Epsilon:new eps, t=32, eps=0.42767785702368233
INFO:ABC:t: 32, eps: 0.42767785702368233.
DEBUG:ABC:Now submitting population 32.
DEBUG:ABC:Population 32 done.
DEBUG:ABC:Total samples up to t = 32: 840214.
INFO:ABC:Acceptance rate: 2000 / 11259 = 1.7764e-01.
DEBUG:Epsilon:new eps, t=33, eps=0.4272103409556444
INFO:ABC:t: 33, eps: 0.4272103409556444.
DEBUG:ABC:Now submitting population 33.
DEBUG:ABC:Population 33 done.
DEBUG:ABC:Total samples up to t = 33: 854209.
INFO:ABC:Acceptance rate: 2000 / 13995 = 1.4291e-01.
DEBUG:Epsilon:new eps, t=34, eps=0.42683015417264375
INFO:ABC:t: 34, eps: 0.42683015417264375.
DEBUG:ABC:Now submitting population 34.
DEBUG:ABC:Population 34 done.
DEBUG:ABC:Total samples up to t = 34: 868660.
INFO:ABC:Acceptance rate: 2000 / 14451 = 1.3840e-01.
DEBUG:Epsilon:new eps, t=35, eps=0.42655086154392674
INFO:ABC:t: 35, eps: 0.42655086154392674.
DEBUG:ABC:Now submitting population 35.
DEBUG:ABC:Population 35 done.
DEBUG:ABC:Total samples up to t = 35: 884119.
INFO:ABC:Acceptance rate: 2000 / 15459 = 1.2937e-01.
DEBUG:Epsilon:new eps, t=36, eps=0.42630783492696955
INFO:ABC:t: 36, eps: 0.42630783492696955.
DEBUG:ABC:Now submitting population 36.
DEBUG:ABC:Population 36 done.
DEBUG:ABC:Total samples up to t = 36: 907482.
INFO:ABC:Acceptance rate: 2000 / 23363 = 8.5605e-02.
DEBUG:Epsilon:new eps, t=37, eps=0.4261289504463212
INFO:ABC:t: 37, eps: 0.4261289504463212.
DEBUG:ABC:Now submitting population 37.
DEBUG:ABC:Population 37 done.
DEBUG:ABC:Total samples up to t = 37: 952664.
INFO:ABC:Acceptance rate: 2000 / 45182 = 4.4265e-02.
DEBUG:Epsilon:new eps, t=38, eps=0.42595630527561346
INFO:ABC:t: 38, eps: 0.42595630527561346.
DEBUG:ABC:Now submitting population 38.
DEBUG:ABC:Population 38 done.
DEBUG:ABC:Total samples up to t = 38: 1005325.
INFO:ABC:Acceptance rate: 2000 / 52661 = 3.7979e-02.
DEBUG:Epsilon:new eps, t=39, eps=0.4258018856386711
INFO:ABC:t: 39, eps: 0.4258018856386711.
DEBUG:ABC:Now submitting population 39.
DEBUG:ABC:Population 39 done.
DEBUG:ABC:Total samples up to t = 39: 1067839.
INFO:ABC:Acceptance rate: 2000 / 62514 = 3.1993e-02.
DEBUG:Epsilon:new eps, t=40, eps=0.42563745668433595
INFO:ABC:t: 40, eps: 0.42563745668433595.
DEBUG:ABC:Now submitting population 40.
DEBUG:ABC:Population 40 done.
DEBUG:ABC:Total samples up to t = 40: 1123261.
INFO:ABC:Acceptance rate: 2000 / 55422 = 3.6087e-02.
DEBUG:Epsilon:new eps, t=41, eps=0.4254995740629113
INFO:ABC:t: 41, eps: 0.4254995740629113.
DEBUG:ABC:Now submitting population 41.
DEBUG:ABC:Population 41 done.
DEBUG:ABC:Total samples up to t = 41: 1157955.
INFO:ABC:Acceptance rate: 2000 / 34694 = 5.7647e-02.
DEBUG:Epsilon:new eps, t=42, eps=0.4253862273938184
INFO:ABC:t: 42, eps: 0.4253862273938184.
DEBUG:ABC:Now submitting population 42.
DEBUG:ABC:Population 42 done.
DEBUG:ABC:Total samples up to t = 42: 1187381.
INFO:ABC:Acceptance rate: 2000 / 29426 = 6.7967e-02.
DEBUG:Epsilon:new eps, t=43, eps=0.42528201968543083
INFO:ABC:t: 43, eps: 0.42528201968543083.
DEBUG:ABC:Now submitting population 43.
DEBUG:ABC:Population 43 done.
DEBUG:ABC:Total samples up to t = 43: 1212725.
INFO:ABC:Acceptance rate: 2000 / 25344 = 7.8914e-02.
DEBUG:Epsilon:new eps, t=44, eps=0.4251961950108309
INFO:ABC:t: 44, eps: 0.4251961950108309.
DEBUG:ABC:Now submitting population 44.
DEBUG:ABC:Population 44 done.
DEBUG:ABC:Total samples up to t = 44: 1240064.
INFO:ABC:Acceptance rate: 2000 / 27339 = 7.3156e-02.
DEBUG:Epsilon:new eps, t=45, eps=0.4251166689549215
INFO:ABC:t: 45, eps: 0.4251166689549215.
DEBUG:ABC:Now submitting population 45.
DEBUG:ABC:Population 45 done.
DEBUG:ABC:Total samples up to t = 45: 1253527.
INFO:ABC:Acceptance rate: 2000 / 13463 = 1.4856e-01.
DEBUG:Epsilon:new eps, t=46, eps=0.4250529684125474
INFO:ABC:t: 46, eps: 0.4250529684125474.
DEBUG:ABC:Now submitting population 46.
DEBUG:ABC:Population 46 done.
DEBUG:ABC:Total samples up to t = 46: 1265876.
INFO:ABC:Acceptance rate: 2000 / 12349 = 1.6196e-01.
DEBUG:Epsilon:new eps, t=47, eps=0.4250012566383052
INFO:ABC:t: 47, eps: 0.4250012566383052.
DEBUG:ABC:Now submitting population 47.
DEBUG:ABC:Population 47 done.
DEBUG:ABC:Total samples up to t = 47: 1278386.
INFO:ABC:Acceptance rate: 2000 / 12510 = 1.5987e-01.
DEBUG:Epsilon:new eps, t=48, eps=0.4249604538479238
INFO:ABC:t: 48, eps: 0.4249604538479238.
DEBUG:ABC:Now submitting population 48.
DEBUG:ABC:Population 48 done.
DEBUG:ABC:Total samples up to t = 48: 1293542.
INFO:ABC:Acceptance rate: 2000 / 15156 = 1.3196e-01.
DEBUG:Epsilon:new eps, t=49, eps=0.42492359977618344
INFO:ABC:t: 49, eps: 0.42492359977618344.
DEBUG:ABC:Now submitting population 49.
DEBUG:ABC:Population 49 done.
DEBUG:ABC:Total samples up to t = 49: 1310714.
INFO:ABC:Acceptance rate: 2000 / 17172 = 1.1647e-01.
DEBUG:Epsilon:new eps, t=50, eps=0.4248875005270672
INFO:ABC:t: 50, eps: 0.4248875005270672.
DEBUG:ABC:Now submitting population 50.
DEBUG:ABC:Population 50 done.
DEBUG:ABC:Total samples up to t = 50: 1338094.
INFO:ABC:Acceptance rate: 2000 / 27380 = 7.3046e-02.
DEBUG:Epsilon:new eps, t=51, eps=0.4248554760452807
INFO:ABC:t: 51, eps: 0.4248554760452807.
DEBUG:ABC:Now submitting population 51.
DEBUG:ABC:Population 51 done.
DEBUG:ABC:Total samples up to t = 51: 1369799.
INFO:ABC:Acceptance rate: 2000 / 31705 = 6.3082e-02.
DEBUG:Epsilon:new eps, t=52, eps=0.4248234624424054
INFO:ABC:t: 52, eps: 0.4248234624424054.
DEBUG:ABC:Now submitting population 52.
DEBUG:ABC:Population 52 done.
DEBUG:ABC:Total samples up to t = 52: 1415295.
INFO:ABC:Acceptance rate: 2000 / 45496 = 4.3960e-02.
DEBUG:Epsilon:new eps, t=53, eps=0.42479477279760786
INFO:ABC:t: 53, eps: 0.42479477279760786.
DEBUG:ABC:Now submitting population 53.
DEBUG:ABC:Population 53 done.
DEBUG:ABC:Total samples up to t = 53: 1489604.
INFO:ABC:Acceptance rate: 2000 / 74309 = 2.6915e-02.
DEBUG:Epsilon:new eps, t=54, eps=0.42476885188604646
INFO:ABC:t: 54, eps: 0.42476885188604646.
DEBUG:ABC:Now submitting population 54.
DEBUG:ABC:Population 54 done.
DEBUG:ABC:Total samples up to t = 54: 1596015.
INFO:ABC:Acceptance rate: 2000 / 106411 = 1.8795e-02.
DEBUG:Epsilon:new eps, t=55, eps=0.4247414185381904
INFO:ABC:t: 55, eps: 0.4247414185381904.
DEBUG:ABC:Now submitting population 55.
DEBUG:ABC:Population 55 done.
DEBUG:ABC:Total samples up to t = 55: 1749120.
INFO:ABC:Acceptance rate: 2000 / 153105 = 1.3063e-02.
DEBUG:Epsilon:new eps, t=56, eps=0.42471677797022983
INFO:ABC:t: 56, eps: 0.42471677797022983.
DEBUG:ABC:Now submitting population 56.
DEBUG:ABC:Population 56 done.
DEBUG:ABC:Total samples up to t = 56: 1969344.
INFO:ABC:Acceptance rate: 2000 / 220224 = 9.0817e-03.
DEBUG:Epsilon:new eps, t=57, eps=0.42469372490894597
INFO:History:Done <ABCSMC(id=1, start_time=2019-12-19 11:31:55.095938, end_time=2019-12-19 19:05:03.510737)>

Analysis of results


In [34]:
history = History(db_path)

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

In [36]:
df.describe()


Out[36]:
name ina.a1_m ina.a2_m ina.a3_m ina.a4_m ina.b1_m ina.b2_m
count 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000
mean -99.875803 0.145998 0.215455 4.868289 1.085277 20.074018
std 0.096275 0.000124 0.001506 2.723653 0.005912 0.069211
min -99.999989 0.145752 0.210474 0.003857 1.067430 19.887022
25% -99.947600 0.145909 0.214507 2.546949 1.081246 20.023639
50% -99.897704 0.145973 0.215409 4.733125 1.084927 20.069865
75% -99.827944 0.146059 0.216456 7.126737 1.089592 20.124146
max -99.359168 0.146663 0.220172 9.999977 1.103101 20.273091

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

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

g = plot_sim_results(modelfile,
                     sakakibara_act,
                     schneider_taum_cou_adjust,
                     df=df, w=w)

plt.tight_layout()



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

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


Fast inactivation gate ($h$) calibration


In [6]:
observations, model, summary_statistics = setup(modelfile,
                                                sakakibara_inact_cou_adjust,
                                                sakakibara_inact_kin_fast_cou_adjust,
                                                sakakibara_rec_fast_cou_adjust)

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

In [8]:
g = plot_sim_results(modelfile,
                     sakakibara_inact_cou_adjust,
                     sakakibara_inact_kin_fast_cou_adjust,
                     sakakibara_rec_fast_cou_adjust)



In [9]:
limits = {'ina.c1_h': (-100, 0),
          'log_ina.a1_h': (-2, 1),
          'ina.a2_h': (0, 50),
          'ina.a3_h': (0, 200),
          'ina.b2_h': (0, 100),
          'ina.b3_h': (0, 50),
          'log_ina.b4_h': (-1, 2),
          'log_ina.b5_h': (-3, 0),
          'log_ina.b6_h': (3, 6),
          'log_ina.b7_h': (-2, 1)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

In [11]:
summary_statistics(model(prior.rvs()))


Out[11]:
{'0': 0.9975858422617265,
 '1': 0.9974627058204689,
 '2': 0.9974870984338249,
 '3': 0.9944592195234383,
 '4': 0.9661312608965233,
 '5': 0.6829092262806274,
 '6': 0.1662969891163373,
 '7': 0.01663217010480187,
 '8': 1.544087931460735e-10,
 '9': 3.8241061891543885e-06,
 '10': -1.3380807686156046e-10,
 '11': 5.017197320804844,
 '12': 3.5947181581467964,
 '13': 3.0629702453271057,
 '14': 2.868538745242184,
 '15': 0.14631954846222775,
 '16': 0.48295832090343416,
 '17': 0.9113064992144112,
 '18': 2.236721091521449,
 '19': 6.323408613241784}

Run ABC calibration


In [12]:
db_path = ("sqlite:///" + os.path.join(tempfile.gettempdir(), "courtemanche_ina_hgate_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 1024 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(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': 0.9298311489289642, '1': 0.9298311489289642, '2': 0.9298311489289642, '3': 0.9298311489289642, '4': 0.9298311489289642, '5': 0.34877521029031383, '6': 0.21674341481487852, '7': 0.24006575158114657, '8': 0.5989135286482826, '9': 0.9298311489289642, '10': 0.9298311489289642, '11': 0.6679914857297752, '12': 1.3151082375304948, '13': 1.8297158087380803, '14': 2.5570356595546513, '15': 2.0456285276437214, '16': 1.901285628941935, '17': 0.7393593039114589, '18': 0.44622757416898995, '19': 0.5843318259435246}
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-12-24 18:12:25.098801, 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 [19]:
history = History(db_path)

In [20]:
history.all_runs()


Out[20]:
[<ABCSMC(id=1, start_time=2019-12-24 18:12:25.098801, end_time=2019-12-25 04:46:04.211283)>]

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

In [22]:
df.describe()


Out[22]:
name ina.a2_h ina.a3_h ina.b2_h ina.b3_h ina.c1_h log_ina.a1_h log_ina.b4_h log_ina.b5_h log_ina.b6_h log_ina.b7_h
count 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000 2000.000000
mean 26.964270 88.508097 24.368373 18.548531 -32.363854 -0.348160 0.307530 -1.462810 4.434213 0.300084
std 5.808081 42.464716 23.564982 10.198298 22.027750 0.688205 0.616413 0.442898 0.810173 0.426152
min 11.383654 0.339995 0.026143 0.153184 -79.268905 -1.991852 -0.999735 -2.992745 3.003545 -0.835733
25% 22.815784 55.654679 8.732021 11.822355 -51.855608 -0.876120 -0.217853 -1.593039 3.733385 -0.001593
50% 26.807161 89.857272 15.919898 15.318698 -28.221729 -0.317193 0.541022 -1.315266 4.407137 0.355084
75% 30.919111 120.014206 29.601389 21.961660 -11.756408 0.182453 0.718996 -1.195146 5.116051 0.639597
max 45.599301 195.259843 99.707189 49.833306 -0.976855 0.992533 1.949545 -0.110172 5.998702 0.998932

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

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

g = plot_sim_results(modelfile,
                     sakakibara_inact_cou_adjust,
                     sakakibara_inact_kin_fast_cou_adjust,
                     sakakibara_rec_fast_cou_adjust,
                     df=df, w=w)

plt.tight_layout()



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

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


Slow inactivation gate ($j$) calibration


In [34]:
observations, model, summary_statistics = setup(modelfile,
                                                sakakibara_inact,
                                                sakakibara_inact_kin_slow,
                                                sakakibara_rec_slow)

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

In [36]:
g = plot_sim_results(modelfile,
                     sakakibara_inact,
                     sakakibara_inact_kin_slow,
                     sakakibara_rec_slow)



In [57]:
limits = {'ina.c1_j': (-100, 0),
          'log_ina.a1_j': (3, 7),
          'log_ina.a2_j': (-2, 2),
          'log_ina.a3_j': (-5, -1),
          'log_ina.a4_j': (-4, 0),
          #'ina.a5_j': (0, 100),
          'ina.a6_j': (0, 1),
          'ina.a7_j': (0, 100),
          #'ina.b1_j': (0, 1.0),
          #'log_ina.b2_j': (-8, -4), # set to zero as original effectively zero
          'ina.b3_j': (0, 1),
          'ina.b4_j': (0, 100),
          'ina.b5_j': (0.0, 1.0),
          'log_ina.b6_j': (-4, 0),
          'ina.b7_j': (0, 1),
          'ina.b8_j': (0, 100)}
prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})

Run ABC calibration


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

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

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


Theoretical minimum population size is 8192 particles

In [61]:
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(8500),
             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.9818869939747233, '1': 0.9818869939747233, '2': 0.9818869939747233, '3': 0.9818869939747233, '4': 0.9818869939747233, '5': 0.3683011084316996, '6': 0.2288776196425918, '7': 0.25350563857514163, '8': 0.6324432182903594, '9': 0.9818869939747233, '10': 0.9818869939747233, '11': 2.700189233430489, '12': 0.4957538129507886, '13': 1.412898366909747, '14': 2.700189233430489, '15': 2.160151386744391, '16': 0.7192132991083816, '17': 0.7852456348032116, '18': 0.3846095367327172, '19': 0.2854129531269287}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=100, quantile_multiplier=1

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

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


INFO:History:Start <ABCSMC(id=9, start_time=2020-01-13 17:12:25.673322, 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 samples up to t = 0: 47298.
INFO:ABC:Acceptance rate: 8500 / 47298 = 1.7971e-01.
DEBUG:Epsilon:new eps, t=1, eps=2.8096148971541752
INFO:ABC:t: 1, eps: 2.8096148971541752.
DEBUG:ABC:Now submitting population 1.

Analysis of results


In [65]:
history = History(db_path)

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

In [67]:
df.describe()


Out[67]:
name ina.a6_j ina.a7_j ina.b3_j ina.b4_j ina.b5_j ina.b7_j ina.b8_j ina.c1_j log_ina.a1_j log_ina.a2_j log_ina.a3_j log_ina.a4_j log_ina.b6_j
count 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000 8500.000000
mean 0.472623 37.068816 0.480921 59.463275 0.389263 0.388710 73.374842 -28.552460 4.949727 0.751182 -3.473260 -2.633527 -2.873703
std 0.271778 24.076817 0.268026 24.815499 0.209265 0.267537 17.645426 19.328287 1.066085 0.658238 0.417799 0.660761 0.633958
min 0.000159 0.002710 0.000063 0.006223 0.002487 0.000146 0.168798 -89.969625 3.000577 -0.820180 -4.993633 -3.999895 -3.998268
25% 0.248080 18.641574 0.258497 41.766668 0.244088 0.154249 64.672158 -37.676529 4.059308 0.233040 -3.698849 -3.150395 -3.383767
50% 0.462635 33.538725 0.478728 62.407472 0.375758 0.339317 77.288463 -25.270756 4.947009 0.759506 -3.412920 -2.562463 -2.914418
75% 0.696230 49.389839 0.694397 79.565410 0.524169 0.589420 86.060323 -14.414774 5.826013 1.275882 -3.210055 -2.080738 -2.415750
max 0.999499 99.915503 0.999966 99.972804 0.999865 0.999914 99.965001 -0.003143 6.998912 1.998630 -1.964447 -1.319365 -0.661527

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

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

g = plot_sim_results(modelfile,
                     sakakibara_inact,
                     sakakibara_inact_kin_slow,
                     sakakibara_rec_slow,
                     df=df, w=w)

plt.tight_layout()



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

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

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



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

In [74]:
from ionchannelABC.visualization import plot_variables

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

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

v = np.arange(-140, 50, 5)

cou_par_map = {'mi': 'ina.m_inf',
               'mt': 'ina.tau_m',
               'hi': 'ina.h_inf',
               'ht': 'ina.tau_h',
               'ji': 'ina.j_inf',
               'jt': 'ina.tau_j'}

f, ax = plot_variables(v, cou_par_map,
                       ["models/courtemanche_ina.mmt"],
                       [cou_par_samples_j],
                       original=True,
                       figshape=(2,3))

plt.tight_layout()


/scratch/cph211/miniconda3/envs/ionchannelABC/lib/python3.7/site-packages/pandas/core/frame.py:7138: 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 [76]:
sns.set_context('poster')

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

v = np.arange(-80, 10, 1)

cou_par_map = {'alpha_j': 'ina.alpha_j',
               'beta_j': 'ina.beta_j'}

f, ax = plot_variables(v, cou_par_map,
                       ["models/courtemanche_ina.mmt"],
                       [cou_par_samples_j],
                       original=True,
                       figshape=(1,2))
#ax[1].set_ylim([0, 0.4])

plt.tight_layout()



In [233]:
cou_par_samples_j[0]


Out[233]:
{'ina.a6_j': 0.41301148704886725,
 'ina.a7_j': 52.6161106141938,
 'ina.b1_j': 0.7831193340019316,
 'ina.b3_j': 0.701468612495651,
 'ina.b4_j': 75.57330479413949,
 'ina.b7_j': 0.8075039406278457,
 'ina.b8_j': 93.04635174044711,
 'ina.c1_j': -7.373498107964862,
 'log_ina.a1_j': 5.996046693845665,
 'log_ina.a2_j': -0.5927309272921197,
 'log_ina.a3_j': -3.0840451102747295,
 'log_ina.a4_j': -2.2256390113967024,
 'log_ina.b2_j': -6.545348838091387,
 'log_ina.b6_j': -2.833325336109481}

In [234]:
def beta_part(V, b1, b2, b3, b4):
    return b1*np.exp(-b2*V)/(1+np.exp(-b3*(V+b4)))

In [235]:
def calc_b1(c1, b2, b3, b4, b5, b6, b7, b8):
    return b5*np.exp(-b6*c1)/(1+np.exp(-b7*(c1+b8)))*(1+np.exp(-b3*(c1+b4)))/np.exp(-b2*c1)

In [236]:
calc_b1(-4.89310260541383, 
        10**-5.850309337197719,
        0.10247765803206027,
        88.10017627036915,
        10**-1.360687617227965,
        10**-0.3200446909034119,
        0.8061796342581548,
        30.695088163032647)


Out[236]:
0.4533174106072827

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

In [268]:
V = np.arange(-100, 50, 1)

In [269]:
beta1 = [beta_part(vi, 
          m.get('ina.b5_j').value(),
          m.get('ina.b6_j').value(),
          m.get('ina.b7_j').value(),
          m.get('ina.b8_j').value()) for vi in V]

In [270]:
beta2 = [beta_part(vi, 
          m.get('ina.b1_j').value(),
          m.get('ina.b2_j').value(),
          m.get('ina.b3_j').value(),
          m.get('ina.b4_j').value()) for vi in V]

In [271]:
plt.plot(V, beta1)
plt.plot(V,beta2)


Out[271]:
[<matplotlib.lines.Line2D at 0x7f54fe9e2210>]

In [ ]:


In [ ]: