ASIC.py

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
import pickle as pk
import os


"""
This file simulates our models of ASIC as isolated channels (i.e. not embedded into a cell, as is done with the NEURON 
simulator).
This allows to explore the basic properties of the channel models, such as their pH-dependent inactivation, recovery 
from inactivation, etc.
The Supplementary figure 6 is mostly generated here.
"""


class ASICType1Chosen:
    """
    Homomeric channel model, based on the data from Baron et al., 2008.
    """
    type = "type1"

    def __init__(self, type="type1", pH05_m=6.46, n_m=1.5, pH05_h=7.3, n_h=4.6, state0=None, basepH=7.4, gmax=1, E_rev=50, Vm=-75):
        self.type = type
        self.pH05_m = pH05_m
        self.n_m = n_m
        self.pH05_h = pH05_h
        self.n_h = n_h
        self.a1_tau_m = 0.3487
        self.a2_tau_m = 234.3
        self.b1_tau_m = 10.14
        self.b2_tau_m = -5.553
        self.c1_tau_m = 6.392
        self.c2_tau_m = 8.595

        self.gamma = 0

        if not state0:
            state0 = [self.m_inf(basepH), self.h_inf(basepH)]
        self.state = state0
        self.gmax = gmax
        self.E_rev = E_rev
        self.Vm = Vm

    def m_inf(self, pH):
        return 1. / (1. + 10**(self.n_m*(pH - self.pH05_m)))

    def h_inf(self, pH):
        return 1.3 / (1 + 10**(-self.n_h*(pH - self.pH05_h)))

    def tau_h(self, pH):
        return 49.196 * np.exp(- 34.682 * (pH - 7.144) ** 2) + (pH - 5) * (4.78 - 0.98) / (9 - 5) + 0.98

    def tau_m(self, pH):
        """Same activation time constant as in Alijevic et al., 2020."""
        return  1 / ( self.a1_tau_m / (1 + np.exp(self.b1_tau_m * (pH - self.c1_tau_m))) + self.a2_tau_m / (1 + np.exp(self.b2_tau_m * (pH - self.c2_tau_m))) ) / 1000

    def set_pHstim(self, fun):
        self.pHstim = fun

    def der(self, state, time, pH=None):
        if type (time) != float:
            state, time = time, state
        m, h = state
        if pH is None:
            pH = self.pHstim(time)
        return np.array([(self.m_inf(pH) - m) / self.tau_m(pH),
                         (self.h_inf(pH) - h) / self.tau_h(pH)])

    def i(self, a, s):
        return self.gmax * a * s * (self.Vm - self.E_rev)

    def grel(self, a, s):
        return a * s

    def run_sim(self, times, out="i"):
        states = integrate(self.der, self.state, times)
        self.state = states[-1]
        if out == "i":
            i = [self.i(*state) for state in states]
        elif out == "grel" or out == "g_rel":
            i = [self.grel(*state) for state in states]
        else:
            raise ValueError
        return np.array(i)


class ASICType2Chosen:
    """
    Heteromeric channel model, based on the data from Baron et al., 2008.
    """
    type = "type2"
    def __init__(self, type="type2", pH05_m=6.03, n_m=1.94, pH05_h=6.74, n_h=3.82, state0=None, basepH=7.4, gmax=1,
                 E_rev=50, Vm=-75):
        self.type = type
        self.pH05_m = pH05_m
        self.n_m = n_m
        self.pH05_h = pH05_h
        self.n_h = n_h
        self.a1_tau_m = 0.3487
        self.a2_tau_m = 234.3
        self.b1_tau_m = 10.14
        self.b2_tau_m = -5.553
        self.c1_tau_m = 6.392
        self.c2_tau_m = 8.595

        self.gamma = 0

        if not state0:
            state0 = [self.m_inf(basepH), self.h_inf(basepH)]
        self.state = state0
        self.gmax = gmax
        self.E_rev = E_rev
        self.Vm = Vm

    def m_inf(self, pH):
        return 1. / (1. + 10**(self.n_m*(pH - self.pH05_m)))

    def h_inf(self, pH):
        return 1 / (1 + 10**(-self.n_h*(pH - self.pH05_h)))

    def tau_h(self, pH):
        a = 42.862
        b = 5.375
        c = 6.600
        d = 1.645
        return a * np.exp(- b * np.abs(pH - c)) + d

    def tau_m(self, pH):
        """Same activation time constant as in Alijevic et al., 2020."""
        return  1 / ( self.a1_tau_m / (1 + np.exp(self.b1_tau_m * (pH - self.c1_tau_m))) + self.a2_tau_m / (1 + np.exp(self.b2_tau_m * (pH - self.c2_tau_m))) ) / 1000

    def set_pHstim(self, fun):
        self.pHstim = fun

    def der(self, state, time, pH=None):
        if type (time) != float:
            state, time = time, state
        m, h = state
        if pH is None:
            pH = self.pHstim(time)
        return np.array([(self.m_inf(pH) - m) / self.tau_m(pH),
                         (self.h_inf(pH) - h) / self.tau_h(pH)])

    def i(self, m, h):
        return self.gmax * m * h * (self.Vm - self.E_rev)

    def grel(self, a, s):
        return a * s

    def run_sim(self, times, out="i"):
        states = integrate(self.der, self.state, times)
        self.state = states[-1]
        if out == "i":
            i = [self.i(*state) for state in states]
        elif out == "grel" or out == "g_rel":
            i = [self.grel(*state) for state in states]
        else:
            raise ValueError
        return np.array(i)


def load_channel(which=1):
    if which == 1:
        chan = ASICType1Chosen()
    elif which == 2:
        chan = ASICType2Chosen()
    return chan


def integrate(ode_fun, state0, time):
    out = solve_ivp(ode_fun, (time[0], time[-1]), state0, t_eval=time, max_step=1) # mxstep=100
    state = np.transpose(out["y"])
    return state


def general_pH(tags):
    """
    tags is a list of the form [(pH1, pH1_duration), (pH2, pH2_duration), ...]
    Generates a pH protocol in which pH1 is applied during pH1_duration, then pH2 during pH2_duration...
    Returns a function that, given a time t, returns the pH for this time in the protocol.
    """
    def fun(t):
        time = 0
        for pH, duration in tags:
            time += duration
            if t <= time:
                return pH
        return pH
    return fun


def get_figS6_protocols():
    return [reactivation_protocol(), inactivation_protocol(), single_drop_protocol()]


def reactivation_protocol():
    """
    Generates the recovery from inactivation protocol of Supplementary figure 6A.
    """
    reac_times = [.5, 1, 2, 4, 8]
    return ("Recovery from inactivation", general_pH,
            [[(7.4, 5), (5, 10), (7.4, t_reac), (5, 15-t_reac)] for t_reac in reac_times],
            ["{}".format(t_reac) for t_reac in reac_times])


def inactivation_protocol():
    """
    Generates the pH-dependent inactivation protocol of Supplementary figure 6B.
    """
    testpHs = [8, 7.6, 7.4, 7, 6.6, 6]
    return ("Activation/Inactivation/SSD", general_pH,
            [[(testpH, 40), (5, 5), (8, 20)] for testpH in testpHs],
            ["{}".format(testpH) for testpH in testpHs])


def single_drop_protocol():
    """
    Generates the pH-dependent inactivation protocol of Supplementary figure 6F.
    """
    testpHs = [6, 7.32]
    return ("Simple drop", general_pH,
            [[(7.4, 40), (testpH, 5), (7.4, 40)] for testpH in testpHs],
            ["{}".format(testpH) for testpH in testpHs])


def current_graphs(funs, folder):
    """
    funs must be a list.
    Each element of the list is a tuple that describes a pH protocol. (Each protocol will be shown in a different figure.)
    The tuple is of the form (protocol_name, genfun, params, labels) where:
        - protocol_name is the name/description of the protocol;
        - gen_fun is a function that, given a parameter set, returns a function t -> pH;
        - params is a list of the different parameter sets for this protocol, to be given to gen_fun (e.g. one parameter
            set describes the protocol for one reactivation time, for the 'reactivation' protocol);
        - labels is a list of labels corresponding to each parameter set.
    The pH protocol will be applied to both type1 and type2 models.
    Because the integration can be slow, the results of the simulation will be saved to pickle files, and the data will
    be loaded instead of simulated if the simulation data is already existing.
    folder is the folder in which the simulated data is to be saved.
    """
    os.makedirs(folder, exist_ok=True)
    for k, (title, genfun, params, labels) in enumerate(funs):
        name = f"rev2_try3_{title.replace(' ', '_').replace('/', '-')}"
        time = np.arange(0, sum(x[1] for x in params[0]), 1e-4)
        fig, axs = plt.subplots(3, 1, figsize=(6, 8), tight_layout=True, gridspec_kw={'height_ratios': [1, 2.5, 2.5]})
        for j, (param, label) in enumerate(zip(params, labels)):
            axs[0].plot(time, np.array([genfun(param)(t) for t in time])-0.015*j, label=label)

            for l, asic_class in enumerate([ASICType1Chosen, ASICType2Chosen]):
                filename = os.path.join(folder, name+f"_{asic_class}_{j}.pickle")
                if os.path.exists(filename):
                    print("Loading")
                    with open(filename, "rb") as f:
                        i = pk.load(f)
                else:
                    print("Simulating")
                    asic = asic_class()
                    asic.set_pHstim(genfun(param))
                    i = asic.run_sim(time, out="grel")
                    with open(filename, "wb") as f:
                        pk.dump(i, f)
                axs[l+1].plot(time, i, label=label)
        axs[0].set_ylabel("Applied pH")
        for ax in axs:
            ax.set_xlabel("Time (s)")
        axs[1].set_ylabel("Homomeric relative conductance m * h")
        axs[2].set_ylabel("Heteromeric relative conductance m * h")
    plt.show()


if __name__ == "__main__":
    current_graphs(get_figS6_protocols(), folder="figS6")