In [1]:

    

# PyABC imports
from pyabc import (ABCSMC, Distribution, RV,
                   History, MedianEpsilon)
from pyabc.populationstrategy import ConstantPopulationSize, AdaptivePopulationSize
from pyabc.epsilon import MedianEpsilon
from pyabc.sampler import MulticoreEvalParallelSampler


    

In [2]:

    

# Custom imports
from ionchannelABC import (ion_channel_sum_stats_calculator,
                           IonChannelAcceptor,
                           IonChannelDistance,
                           EfficientMultivariateNormalTransition,
                           calculate_parameter_sensitivity,
                           plot_parameter_sensitivity,
                           plot_regression_fit,
                           plot_parameters_kde,
                           plot_sim_results,
                           plot_distance_weights)


    

    




INFO:myokit:Loading Myokit version 1.27.4

In [3]:

    

# Other necessary imports
import numpy as np
import subprocess
import pandas as pd
import io
import os
import tempfile


    

In [4]:

    

# Plotting imports
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
%config Inline.Backend.figure_format = 'retina'


    

In [5]:

    

from channels.ito_generic import ito as model
model.sample({})


    

    
Out[5]:







  

    

      

      
x
      
y
      
exp
    
  
  

    

      
0
      
-40
      
1.134054
      
0
    
    

      
1
      
-30
      
2.073325
      
0
    
    

      
2
      
-20
      
3.202193
      
0
    
    

      
3
      
-10
      
3.338438
      
0
    
    

      
4
      
0
      
5.131921
      
0
    
    

      
5
      
10
      
5.945935
      
0
    
    

      
6
      
20
      
6.609532
      
0
    
    

      
7
      
30
      
7.233689
      
0
    
    

      
8
      
40
      
7.839728
      
0
    
    

      
9
      
50
      
8.437684
      
0
    
    

      
10
      
60
      
9.032094
      
0
    
    

      
0
      
0
      
0.855227
      
1
    
    

      
1
      
10
      
0.840949
      
1
    
    

      
2
      
20
      
0.818564
      
1
    
    

      
3
      
30
      
0.825545
      
1
    
    

      
4
      
40
      
0.828194
      
1
    
    

      
5
      
50
      
0.829217
      
1
    
    

      
6
      
60
      
0.829616
      
1
    
    

      
0
      
-100
      
1.000000
      
2
    
    

      
1
      
-90
      
0.992242
      
2
    
    

      
2
      
-80
      
0.974369
      
2
    
    

      
3
      
-70
      
0.935473
      
2
    
    

      
4
      
-60
      
0.852997
      
2
    
    

      
5
      
-50
      
0.706863
      
2
    
    

      
6
      
-40
      
0.503757
      
2
    
    

      
7
      
-30
      
0.302487
      
2
    
    

      
8
      
-20
      
0.158172
      
2
    
    

      
9
      
-10
      
0.076150
      
2
    
  

In [6]:

    

measurements = model.get_experiment_data()
obs = measurements.to_dict()['y']
exp = measurements.to_dict()['exp']
errs = measurements.to_dict()['errs']


    

In [7]:

    

limits = dict(g_to=(0, 1),
              Vhalf_x=(-100,100),
              k_x=(0,100),
              c_bx=(0,10),
              c_ax=(0,100),
              sigma_x=(0,50),
              Vmax_x=(-100,100),
              Vhalf_y=(-100,100),
              k_y=(-100,0),
              c_by=(0,100),
              c_ay=(0,100),
              sigma_y=(0,100),
              Vmax_y=(-100,100))


    

In [8]:

    

parameters = ['ito.'+k for k in limits.keys()]


    

In [9]:

    

distance_fn=IonChannelDistance(
    obs=obs,
    exp_map=exp,
    err_bars=errs,
    err_th=0.1)


    

In [10]:

    

sns.set_context('talk')
g = plot_distance_weights(model, distance_fn)


    

In [11]:

    

g.savefig('results/ito-generic/dist-weights.pdf')


    

In [15]:

    

fitted, regression_fit, r2 = calculate_parameter_sensitivity(
    model,
    parameters,
    distance_fn,
    sigma=0.1,
    n_samples=1000)


    

In [16]:

    

grid1 = plot_parameter_sensitivity(fitted, plot_cutoff=0.05)


    

In [17]:

    

grid2 = plot_regression_fit(regression_fit, r2)


    

In [18]:

    

grid1.savefig('results/ito-generic/sensitivity.pdf')
grid2.savefig('results/ito-generic/sensitivity_fit.pdf')


    

In [19]:

    

prior = Distribution(**{key: RV("uniform", a, b - a)
                        for key, (a,b) in limits.items()})


    

In [20]:

    

db_path = ('sqlite:///' + 
           os.path.join(tempfile.gettempdir(), "hl-1_ito-generic.db"))
print(db_path)


    

    




sqlite:////scratch/cph211/tmp/hl-1_ito-generic.db

In [21]:

    

# Let's log all the sh!t
import logging
logging.basicConfig()
abc_logger = logging.getLogger('ABC')
abc_logger.setLevel(logging.DEBUG)
eps_logger = logging.getLogger('Epsilon')
eps_logger.setLevel(logging.DEBUG)
cv_logger = logging.getLogger('CV Estimation')
cv_logger.setLevel(logging.DEBUG)


    

In [22]:

    

abc = ABCSMC(models=model,
             parameter_priors=prior,
             distance_function=IonChannelDistance(
                 obs=obs,
                 exp_map=exp,
                 err_bars=errs,
                 err_th=0.1),
             population_size=ConstantPopulationSize(2000),
             #population_size=AdaptivePopulationSize(
             #    start_nr_particles=2500,
             #    mean_cv=0.5,
             #    max_population_size=5000,
             #    min_population_size=500),
             summary_statistics=ion_channel_sum_stats_calculator,
             transitions=EfficientMultivariateNormalTransition(),
             eps=MedianEpsilon(),
             sampler=MulticoreEvalParallelSampler(n_procs=6),
             acceptor=IonChannelAcceptor())


    

    




DEBUG:ABC:ion channel weights: {0: 0.28000483001869525, 1: 0.28000483001869525, 2: 0.28000483001869525, 3: 0.28000483001869525, 4: 0.28000483001869525, 5: 0.28000483001869525, 6: 0.28000483001869525, 7: 0.25303936896267654, 8: 0.22088811308980946, 9: 0.17954942768846782, 10: 0.15514335966426648, 11: 0.24096649623380578, 12: 0.34990528496385936, 13: 0.32069909158384313, 14: 0.38470941572561096, 15: 0.4046907051929525, 16: 0.42717066503692636, 17: 0.38491728580794116, 18: 3.710177692724231, 19: 0.6379695788708737, 20: 0.8437662172163188, 21: 2.906305859300656, 22: 2.377886612155085, 23: 2.1797293944754856, 24: 3.710177692724231, 25: 1.4531529296503196, 26: 1.188943306077537, 27: 3.710177692724231}
DEBUG:Epsilon:init quantile_epsilon initial_epsilon=from_sample, quantile_multiplier=1

In [23]:

    

abc_id = abc.new(db_path, obs)


    

    




INFO:History:Start <ABCSMC(id=1, start_time=2018-11-16 17:34:24.800545, end_time=None)>
INFO:Epsilon:initial epsilon is 14.182916273589099

In [24]:

    

history = abc.run(minimum_epsilon=0.01, max_nr_populations=30, min_acceptance_rate=0.01)


    

    




INFO:ABC:t:0 eps:14.182916273589099
DEBUG:ABC:now submitting population 0
DEBUG:ABC:population 0 done
DEBUG:ABC:
total nr simulations up to t =0 is 4038
DEBUG:Epsilon:new eps, t=1, eps=7.003082317131013
INFO:ABC:t:1 eps:7.003082317131013
DEBUG:ABC:now submitting population 1
DEBUG:ABC:population 1 done
DEBUG:ABC:
total nr simulations up to t =1 is 8868
DEBUG:Epsilon:new eps, t=2, eps=5.437646759682292
INFO:ABC:t:2 eps:5.437646759682292
DEBUG:ABC:now submitting population 2
DEBUG:ABC:population 2 done
DEBUG:ABC:
total nr simulations up to t =2 is 14027
DEBUG:Epsilon:new eps, t=3, eps=4.92812636612909
INFO:ABC:t:3 eps:4.92812636612909
DEBUG:ABC:now submitting population 3
DEBUG:ABC:population 3 done
DEBUG:ABC:
total nr simulations up to t =3 is 21092
DEBUG:Epsilon:new eps, t=4, eps=4.4283151717696105
INFO:ABC:t:4 eps:4.4283151717696105
DEBUG:ABC:now submitting population 4
DEBUG:ABC:population 4 done
DEBUG:ABC:
total nr simulations up to t =4 is 29135
DEBUG:Epsilon:new eps, t=5, eps=4.074411861421005
INFO:ABC:t:5 eps:4.074411861421005
DEBUG:ABC:now submitting population 5
DEBUG:ABC:population 5 done
DEBUG:ABC:
total nr simulations up to t =5 is 38815
DEBUG:Epsilon:new eps, t=6, eps=3.802959479143839
INFO:ABC:t:6 eps:3.802959479143839
DEBUG:ABC:now submitting population 6
DEBUG:ABC:population 6 done
DEBUG:ABC:
total nr simulations up to t =6 is 53625
DEBUG:Epsilon:new eps, t=7, eps=3.5953655956510806
INFO:ABC:t:7 eps:3.5953655956510806
DEBUG:ABC:now submitting population 7
DEBUG:ABC:population 7 done
DEBUG:ABC:
total nr simulations up to t =7 is 68595
DEBUG:Epsilon:new eps, t=8, eps=3.3590572270669656
INFO:ABC:t:8 eps:3.3590572270669656
DEBUG:ABC:now submitting population 8
DEBUG:ABC:population 8 done
DEBUG:ABC:
total nr simulations up to t =8 is 81510
DEBUG:Epsilon:new eps, t=9, eps=3.198701603108157
INFO:ABC:t:9 eps:3.198701603108157
DEBUG:ABC:now submitting population 9
DEBUG:ABC:population 9 done
DEBUG:ABC:
total nr simulations up to t =9 is 93333
DEBUG:Epsilon:new eps, t=10, eps=2.9956888524822145
INFO:ABC:t:10 eps:2.9956888524822145
DEBUG:ABC:now submitting population 10
DEBUG:ABC:population 10 done
DEBUG:ABC:
total nr simulations up to t =10 is 106248
DEBUG:Epsilon:new eps, t=11, eps=2.859446330753542
INFO:ABC:t:11 eps:2.859446330753542
DEBUG:ABC:now submitting population 11
DEBUG:ABC:population 11 done
DEBUG:ABC:
total nr simulations up to t =11 is 123116
DEBUG:Epsilon:new eps, t=12, eps=2.7205699469238107
INFO:ABC:t:12 eps:2.7205699469238107
DEBUG:ABC:now submitting population 12
DEBUG:ABC:population 12 done
DEBUG:ABC:
total nr simulations up to t =12 is 142105
DEBUG:Epsilon:new eps, t=13, eps=2.5491609207197055
INFO:ABC:t:13 eps:2.5491609207197055
DEBUG:ABC:now submitting population 13
DEBUG:ABC:population 13 done
DEBUG:ABC:
total nr simulations up to t =13 is 160640
DEBUG:Epsilon:new eps, t=14, eps=2.4117434115817953
INFO:ABC:t:14 eps:2.4117434115817953
DEBUG:ABC:now submitting population 14
DEBUG:ABC:population 14 done
DEBUG:ABC:
total nr simulations up to t =14 is 174550
DEBUG:Epsilon:new eps, t=15, eps=2.240575346314321
INFO:ABC:t:15 eps:2.240575346314321
DEBUG:ABC:now submitting population 15
DEBUG:ABC:population 15 done
DEBUG:ABC:
total nr simulations up to t =15 is 191585
DEBUG:Epsilon:new eps, t=16, eps=2.1116953392020448
INFO:ABC:t:16 eps:2.1116953392020448
DEBUG:ABC:now submitting population 16
DEBUG:ABC:population 16 done
DEBUG:ABC:
total nr simulations up to t =16 is 211252
DEBUG:Epsilon:new eps, t=17, eps=1.9476333926982359
INFO:ABC:t:17 eps:1.9476333926982359
DEBUG:ABC:now submitting population 17
DEBUG:ABC:population 17 done
DEBUG:ABC:
total nr simulations up to t =17 is 230803
DEBUG:Epsilon:new eps, t=18, eps=1.704354956202357
INFO:ABC:t:18 eps:1.704354956202357
DEBUG:ABC:now submitting population 18
DEBUG:ABC:population 18 done
DEBUG:ABC:
total nr simulations up to t =18 is 276820
DEBUG:Epsilon:new eps, t=19, eps=1.548059143680667
INFO:ABC:t:19 eps:1.548059143680667
DEBUG:ABC:now submitting population 19
DEBUG:ABC:population 19 done
DEBUG:ABC:
total nr simulations up to t =19 is 291957
DEBUG:Epsilon:new eps, t=20, eps=1.4364933320731321
INFO:ABC:t:20 eps:1.4364933320731321
DEBUG:ABC:now submitting population 20
DEBUG:ABC:population 20 done
DEBUG:ABC:
total nr simulations up to t =20 is 308625
DEBUG:Epsilon:new eps, t=21, eps=1.3138367625803244
INFO:ABC:t:21 eps:1.3138367625803244
DEBUG:ABC:now submitting population 21
DEBUG:ABC:population 21 done
DEBUG:ABC:
total nr simulations up to t =21 is 324922
DEBUG:Epsilon:new eps, t=22, eps=1.213371803509123
INFO:ABC:t:22 eps:1.213371803509123
DEBUG:ABC:now submitting population 22
DEBUG:ABC:population 22 done
DEBUG:ABC:
total nr simulations up to t =22 is 341112
DEBUG:Epsilon:new eps, t=23, eps=1.1339092975746319
INFO:ABC:t:23 eps:1.1339092975746319
DEBUG:ABC:now submitting population 23
DEBUG:ABC:population 23 done
DEBUG:ABC:
total nr simulations up to t =23 is 358174
DEBUG:Epsilon:new eps, t=24, eps=1.066186356032324
INFO:ABC:t:24 eps:1.066186356032324
DEBUG:ABC:now submitting population 24
DEBUG:ABC:population 24 done
DEBUG:ABC:
total nr simulations up to t =24 is 380812
DEBUG:Epsilon:new eps, t=25, eps=0.9822637121279986
INFO:ABC:t:25 eps:0.9822637121279986
DEBUG:ABC:now submitting population 25
DEBUG:ABC:population 25 done
DEBUG:ABC:
total nr simulations up to t =25 is 405074
DEBUG:Epsilon:new eps, t=26, eps=0.9092689664418522
INFO:ABC:t:26 eps:0.9092689664418522
DEBUG:ABC:now submitting population 26
DEBUG:ABC:population 26 done
DEBUG:ABC:
total nr simulations up to t =26 is 430445
DEBUG:Epsilon:new eps, t=27, eps=0.8483403493276322
INFO:ABC:t:27 eps:0.8483403493276322
DEBUG:ABC:now submitting population 27
DEBUG:ABC:population 27 done
DEBUG:ABC:
total nr simulations up to t =27 is 472915
DEBUG:Epsilon:new eps, t=28, eps=0.8013815486857275
INFO:ABC:t:28 eps:0.8013815486857275
DEBUG:ABC:now submitting population 28
DEBUG:ABC:population 28 done
DEBUG:ABC:
total nr simulations up to t =28 is 495001
DEBUG:Epsilon:new eps, t=29, eps=0.7475295926330194
INFO:ABC:t:29 eps:0.7475295926330194
DEBUG:ABC:now submitting population 29
DEBUG:ABC:population 29 done
DEBUG:ABC:
total nr simulations up to t =29 is 521882
DEBUG:Epsilon:new eps, t=30, eps=0.6976099126606371
INFO:History:Done <ABCSMC(id=1, start_time=2018-11-16 17:34:24.800545, end_time=2018-11-16 21:35:25.297968)>

In [ ]:

    

history = abc.run(minimum_epsilon=0.01, max_nr_populations=30, min_acceptance_rate=0.01)


    

    




INFO:ABC:t:30 eps:0.6976099126606371
DEBUG:ABC:now submitting population 30
DEBUG:ABC:population 30 done
DEBUG:ABC:
total nr simulations up to t =30 is 547948
DEBUG:Epsilon:new eps, t=31, eps=0.667321402455044
INFO:ABC:t:31 eps:0.667321402455044
DEBUG:ABC:now submitting population 31
DEBUG:ABC:population 31 done
DEBUG:ABC:
total nr simulations up to t =31 is 578745
DEBUG:Epsilon:new eps, t=32, eps=0.6354312238781075
INFO:ABC:t:32 eps:0.6354312238781075
DEBUG:ABC:now submitting population 32
DEBUG:ABC:population 32 done
DEBUG:ABC:
total nr simulations up to t =32 is 604634
DEBUG:Epsilon:new eps, t=33, eps=0.5876082489546789
INFO:ABC:t:33 eps:0.5876082489546789
DEBUG:ABC:now submitting population 33
DEBUG:ABC:population 33 done
DEBUG:ABC:
total nr simulations up to t =33 is 641234
DEBUG:Epsilon:new eps, t=34, eps=0.5543424367982401
INFO:ABC:t:34 eps:0.5543424367982401
DEBUG:ABC:now submitting population 34
DEBUG:ABC:population 34 done
DEBUG:ABC:
total nr simulations up to t =34 is 666774
DEBUG:Epsilon:new eps, t=35, eps=0.5248461771113363
INFO:ABC:t:35 eps:0.5248461771113363
DEBUG:ABC:now submitting population 35
DEBUG:ABC:population 35 done
DEBUG:ABC:
total nr simulations up to t =35 is 692448
DEBUG:Epsilon:new eps, t=36, eps=0.49735755508451307
INFO:ABC:t:36 eps:0.49735755508451307
DEBUG:ABC:now submitting population 36
DEBUG:ABC:population 36 done
DEBUG:ABC:
total nr simulations up to t =36 is 715786
DEBUG:Epsilon:new eps, t=37, eps=0.46759068268909076
INFO:ABC:t:37 eps:0.46759068268909076
DEBUG:ABC:now submitting population 37
DEBUG:ABC:population 37 done
DEBUG:ABC:
total nr simulations up to t =37 is 743859
DEBUG:Epsilon:new eps, t=38, eps=0.44447724342925166
INFO:ABC:t:38 eps:0.44447724342925166
DEBUG:ABC:now submitting population 38

In [26]:

    

db_path = 'sqlite:////scratch/cph211/tmp/hl-1_ito-generic.db'
history = History(db_path)
history.all_runs()


    

    
Out[26]:





[<ABCSMC(id=1, start_time=2018-11-16 17:34:24.800545, end_time=2018-11-17 06:49:09.662319)>]

In [27]:

    

history.id = 1


    

In [28]:

    

sns.set_context('talk')


    

In [29]:

    

evolution = history.get_all_populations()
grid = sns.relplot(x='t', y='epsilon', size='samples', data=evolution[evolution.t>=0])
grid.savefig('results/ito-generic/eps_evolution.pdf')


    

In [30]:

    

df, w = history.get_distribution(m=0)


    

In [31]:

    

g = plot_parameters_kde(df, w, limits, aspect=5, height=1.1)


    

In [32]:

    

g.savefig('results/ito-generic/parameters_kde.pdf')


    

In [33]:

    

# Generate parameter samples
n_samples = 100
df, w = history.get_distribution(m=0)
th_samples = df.sample(n=n_samples, weights=w, replace=True).to_dict(orient='records')


    

In [34]:

    

# Generate sim results samples
samples = pd.DataFrame({})
for i, th in enumerate(th_samples):
    output = model.sample(pars=th, n_x=50)
    output['sample'] = i
    output['distribution'] = 'post'
    samples = samples.append(output, ignore_index=True)


    

In [35]:

    

from ionchannelABC import plot_sim_results
sns.set_context('talk')
g = plot_sim_results(samples, obs=measurements)

# Set axis labels
xlabels = ["voltage, mV", "voltage, mV", "voltage, mV"]
ylabels = ["current density, pA/pF", "activation time constant", "inactivation"]
for ax, xl in zip(g.axes.flatten(), xlabels):
    ax.set_xlabel(xl)
for ax, yl in zip(g.axes.flatten(), ylabels):
    ax.set_ylabel(yl)


    

In [36]:

    

g.savefig('results/ito-generic/ito_sim_results.pdf')


    

In [19]:

    

# Fit activation time constants
grouped = samples[samples['exp']==1].groupby('sample')


    

In [20]:

    

from scipy.optimize import curve_fit
def fit_single_exp(group):
    def single_exp(V, a, b, c):
        return a + b*(np.exp(-V/c))
    guess = (1, 10, 30)
    popt, _ = curve_fit(single_exp, group.x, group.y, guess)
    return popt
output = grouped.apply(fit_single_exp).apply(pd.Series)


    

In [21]:

    

print(output.mean())
print(output.std())


    

    




0     0.693516
1    11.346714
2    26.009547
dtype: float64
0    0.067142
1    0.064778
2    0.400484
dtype: float64

In [22]:

    

import scipy.stats as st
a = output[0].tolist()
rv = st.rv_discrete(values=(a, [1/len(a),]*len(a)))
print("median: {}".format(rv.median()))
print("95% CI: {}".format(rv.interval(0.95)))


    

    




median: 0.698940394909814
95% CI: (0.5624895996636707, 0.7990268921811594)

In [23]:

    

b = output[1].tolist()
rv = st.rv_discrete(values=(b, [1/len(b),]*len(b)))
print("median: {}".format(rv.median()))
print("95% CI: {}".format(rv.interval(0.95)))


    

    




median: 11.346714970616826
95% CI: (11.218115990734853, 11.479266371921778)

In [25]:

    

c = output[2].tolist()
rv = st.rv_discrete(values=(c, [1/len(c),]*len(c)))
print("median: {}".format(rv.median()))
print("95% CI: {}".format(rv.interval(0.95)))


    

    




median: 26.00068489921095
95% CI: (25.23582943151088, 26.870679423724155)

In [26]:

    

# Steady-state inactivation fit to Boltzmann
grouped = samples[samples['exp']==2].groupby('sample')
def fit_boltzmann(group):
    def boltzmann(V, Vhalf, Shalf):
        return 1/(1+np.exp((V-Vhalf)/Shalf))
    guess = (50, 10)
    popt, _ = curve_fit(boltzmann, group.x, group.y, guess)
    return popt
output = grouped.apply(fit_boltzmann).apply(pd.Series)


    

In [27]:

    

print(output.mean())
print(output.std())


    

    




0   -24.478147
1     4.032313
dtype: float64
0    0.136882
1    0.079504
dtype: float64

In [28]:

    

Vhalf = output[0].tolist()
rv = st.rv_discrete(values=(Vhalf, [1/len(Vhalf),]*len(Vhalf)))
print("median: {}".format(rv.median()))
print("95% CI: {}".format(rv.interval(0.95)))


    

    




median: -24.48606536212124
95% CI: (-24.767748128562673, -24.23665583915842)

In [30]:

    

slope = output[1].tolist()
rv = st.rv_discrete(values=(slope, [1/len(slope),]*len(slope)))
print("median: {}".format(rv.median()))
print("95% CI: {}".format(rv.interval(0.95)))


    

    




median: 4.041529219873673
95% CI: (3.869016395723247, 4.170651871570859)

In [34]:

    

from ionchannelABC import plot_distance_weights
grid = plot_distance_weights(model, distance_fn)


    

In [35]:

    

grid.savefig('results/ito/dist_weights.pdf')


    

In [ ]: