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
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
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'
Combine model and experiments to produce:
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)
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())))
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))
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())
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)
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)
In [34]:
history = History(db_path)
In [35]:
df, w = history.get_distribution(m=0)
In [36]:
df.describe()
Out[36]:
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()
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]:
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))
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())
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)
In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [19]:
history = History(db_path)
In [20]:
history.all_runs()
Out[20]:
In [21]:
df, w = history.get_distribution(m=0)
In [22]:
df.describe()
Out[22]:
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()
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()})
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))
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())
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)
In [ ]:
history = abc.run(minimum_epsilon=0., max_nr_populations=100, min_acceptance_rate=0.01)
In [65]:
history = History(db_path)
In [66]:
df, w = history.get_distribution(m=0)
In [67]:
df.describe()
Out[67]:
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()
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]:
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]:
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]:
In [ ]:
In [ ]: