Loading WiChR_model.py 0 → 100644 +236 −0 Original line number Diff line number Diff line import numpy as np from scipy import integrate import matplotlib.pyplot as plt def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, flag_ode=0): """ :param t: time in ms :param X: list containing the states in the model (X=[C1, O1, O2, C1p, nai, ki]). Units of nai and ki are M :param v: transmembrane voltage in mV :param WiChR_g: WiChR conductance in nS/pF :param V_cell: cell volume in l :param C_cell: cell capacitance in F :param irr_times: list containing the illumination start times :param opt_dur: duration of the light pulse in ms :param opt_irr_max: irradiance in mW/mm^2 :param model: 'one-branch' or 'two-branch' :param flag_ode: 0 if the state updates should be outputted, 1 if the WiChR current should be outputted :return: states updates [1/ms] or WiChR current [pA/pF] """ # constants R = 8314.0 # mJ/mol/K T = 310.0 # K F = 96485.0 # C/mol if v != 0: FvRT = F * (v) / (R * T) else: FvRT = F * (1e-2) / (R * T) # ion concentrations nao = 110e-3 # M ko = 1e-3 # M # patch pipette parameters D_K = 187 # Diffusion coefficient in 10e-7 cm^2/s D_Na = 131 # Diffusion coefficient in 10e-7 cm^2/s tau_K = 6000 # ms tau_Na = D_K/D_Na * tau_K # ms na_pp = 1e-3 # M k_pp = 110e-3 # M # illumination opt_irr = 0 for opt_start in irr_times: if opt_start <= t < opt_start + opt_dur: opt_irr = opt_irr_max # states nai = X[0] # M ki = X[1] # M WiChR_C1 = X[2] WiChR_C1p = X[3] WiChR_O1 = X[4] WiChR_O2 = X[5] WiChR_O3 = X[6] WiChR_C2 = X[7] WiChR_O4 = X[8] if model == 'one-branch': # WiChR parameters PKPNa1 = 62.5 PKPNa2 = 52.4 WiChR_g2g1 = 0.904 WiChR_k = 0.229 # mm^2/mW/ms WiChR_a = 5.93 # 1/ms WiChR_e = 0.0033 # 1/ms/mV WiChR_d = 0.0035 # 1/ms # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R*T/F * np.log((PKPNa1*ko + nao)/(PKPNa1*ki + nai)) EK2 = R*T/F * np.log((PKPNa2*ko + nao)/(PKPNa2*ki + nai)) # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 WiChR_dC1pO1 = WiChR_a * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v+150) * WiChR_O1 WiChR_dO2C1 = WiChR_d * WiChR_O2 WiChR_dC1p = WiChR_dC1C1p - WiChR_dC1pO1 WiChR_dC1 = WiChR_dO2C1 - WiChR_dC1C1p WiChR_dO1 = WiChR_dC1pO1 - WiChR_dO1O2 WiChR_dO2 = WiChR_dO1O2 - WiChR_dO2C1 WiChR_dO3 = 0 WiChR_dC2 = 0 WiChR_dO4 = 0 # WiChR current I_WiChR1 = (WiChR_g * WiChR_O1) * (v - EK1) # pA/pF I_WiChR2 = (WiChR_g2g1 * WiChR_g * WiChR_O2) * (v - EK2) # pA/pF I_WiChR = I_WiChR1 + I_WiChR2 # pA/pF # Calculation of the K and Na current through WiChR (optional to use) WiChR_PNa1 = (I_WiChR1) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa1 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa2 = (I_WiChR2) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa2 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_I_K = FvRT * F * PKPNa1 * (WiChR_PNa1) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa2 * (WiChR_PNa2) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2) * (nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA elif model == 'two-branch': # WiChR parameters PKPNa1 = 79.1 PKPNa2 = 58.5 PKPNa3 = 54.1 PKPNa4 = 27.7 WiChR_g2g1 = 0.93 WiChR_g3g1 = 1.06 WiChR_g4g1 = 0.89 WiChR_k = 0.33 # mm^2/mW/ms if opt_irr == 0: WiChR_k1 = 0 else: WiChR_k1 = 24.4/(1+(50/opt_irr)**1.2) # 1/ms WiChR_k2 = 0.01 # mm^2/mW/ms WiChR_a = 6.0 # 1/ms WiChR_e = 0.0031 # 1/ms/mV WiChR_e1 = 0.0042 # 1/ms WiChR_d = 0.00315 # 1/ms WiChR_d1 = 0.0001 # 1/ms WiChR_d2 = 0.0044 # 1/ms WiChR_r = 0.00002 # 1/ms # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R * T / F * np.log((PKPNa1 * ko + nao) / (PKPNa1 * ki + nai)) EK2 = R * T / F * np.log((PKPNa2 * ko + nao) / (PKPNa2 * ki + nai)) EK3 = R * T / F * np.log((PKPNa3 * ko + nao) / (PKPNa3 * ki + nai)) EK4 = R * T / F * np.log((PKPNa4 * ko + nao) / (PKPNa4 * ki + nai)) # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 WiChR_dC1pO1 = WiChR_a * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v + 150) * WiChR_O1 WiChR_dO2C1 = WiChR_d * WiChR_O2 WiChR_dO2O3 = WiChR_k1 * WiChR_O2 WiChR_dO3O2 = WiChR_e1 * WiChR_O3 WiChR_dO3C2 = WiChR_d1 * WiChR_O3 WiChR_dC2O4 = WiChR_k2 * opt_irr * WiChR_C2 WiChR_dO4C2 = WiChR_d2 * WiChR_O4 WiChR_dC2C1 = WiChR_r * WiChR_C2 WiChR_dC1 = WiChR_dO2C1 + WiChR_dC2C1 - WiChR_dC1C1p WiChR_dC1p = WiChR_dC1C1p - WiChR_dC1pO1 WiChR_dO1 = WiChR_dC1pO1 - WiChR_dO1O2 WiChR_dO2 = WiChR_dO1O2 + WiChR_dO3O2 - WiChR_dO2C1 - WiChR_dO2O3 WiChR_dO3 = WiChR_dO2O3 - WiChR_dO3O2 - WiChR_dO3C2 WiChR_dC2 = WiChR_dO3C2 + WiChR_dO4C2 - WiChR_dC2O4 - WiChR_dC2C1 WiChR_dO4 = WiChR_dC2O4 - WiChR_dO4C2 # WiChR current I_WiChR1 = (WiChR_g * WiChR_O1) * (v - EK1) # pA/pF I_WiChR2 = (WiChR_g2g1 * WiChR_g * WiChR_O2) * (v - EK2) # pA/pF I_WiChR3 = (WiChR_g3g1 * WiChR_g * WiChR_O3) * (v - EK3) # pA/pF I_WiChR4 = (WiChR_g4g1 * WiChR_g * WiChR_O4) * (v - EK4) # pA/pF I_WiChR = I_WiChR1 + I_WiChR2 + I_WiChR3 + I_WiChR4 # pA/pF # Calculation of the K and Na current through WiChR (optional to use) WiChR_PNa1 = (I_WiChR1) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa1 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa2 = (I_WiChR2) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa2 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa3 = (I_WiChR3) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa3 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa4 = (I_WiChR4) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa4 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_I_K = FvRT * F * PKPNa1 * (WiChR_PNa1) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa2 * (WiChR_PNa2) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa3 * (WiChR_PNa3) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa4 * (WiChR_PNa4) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2 + WiChR_PNa3 + WiChR_PNa4) * (nai - nao * np.exp(-FvRT)) \ / (1 - np.exp(-FvRT)) # pA Q_K = -WiChR_I_K * C_cell # Total charge in pC n_current_K = Q_K / F # Moles of K due to current delta_concentration_current_K = n_current_K / V_cell / 1000 # Concentration change in M/ms Q_Na = -WiChR_I_Na * C_cell # Total charge in pC n_current_Na = Q_Na / F # Moles of Na due to current delta_concentration_current_Na = n_current_Na / V_cell / 1000 # Concentration change in M/ms # Update concentration due to diffusion delta_pipette_K = - 1 / tau_K * (ki - k_pp) # M/ms delta_pipette_Na = - 1 / tau_Na * (nai - na_pp) # M/ms dK = delta_concentration_current_K + delta_pipette_K # M/ms dNa = delta_concentration_current_Na + delta_pipette_Na # M/ms if flag_ode == 0: output = np.array([dNa, dK, WiChR_dC1, WiChR_dC1p, WiChR_dO1, WiChR_dO2, WiChR_dO3, WiChR_dC2, WiChR_dO4]) else: output = I_WiChR return output def simulate_WiChR(v, opt_irr_max, opt_dur, WiChR_g, V_cell, C_cell, model='one-branch', irr_times=[0], dt_sim=1, off_dur=2000, X0=[1e-3,110e-3,1,0,0,0,0,0,0]): """ Simulate the behavior of WiChR. :param v: transmembrane voltage in mV :param opt_irr_max: irradiance in mW/mm^2 :param opt_dur: duration of the light pulse in ms :param WiChR_g: WiChR conductance in nS/pF :param V_cell: cell volume in l :param C_cell: cell capacitance in F :param model: 'one-branch' or 'two-branch' photocycle model :param irr_times: list containing the illumination start times :param dt_sim: time step of the simulation in ms :param off_dur: duration after the light pulse that should be simulated :param X0: initial conditions of the WiChR model states :return: WiChR current in pA/pF """ time = np.arange(0, irr_times[-1] + opt_dur + off_dur, dt_sim) X = integrate.odeint(WiChR_model, X0, time, tfirst=True, args=(v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model), hmax=1) I_WiChR = [] for i in range(time.size): I_WiChR.append( WiChR_model(time[i], X[i, :], v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, flag_ode=1)) return time, np.array(I_WiChR) Loading
WiChR_model.py 0 → 100644 +236 −0 Original line number Diff line number Diff line import numpy as np from scipy import integrate import matplotlib.pyplot as plt def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, flag_ode=0): """ :param t: time in ms :param X: list containing the states in the model (X=[C1, O1, O2, C1p, nai, ki]). Units of nai and ki are M :param v: transmembrane voltage in mV :param WiChR_g: WiChR conductance in nS/pF :param V_cell: cell volume in l :param C_cell: cell capacitance in F :param irr_times: list containing the illumination start times :param opt_dur: duration of the light pulse in ms :param opt_irr_max: irradiance in mW/mm^2 :param model: 'one-branch' or 'two-branch' :param flag_ode: 0 if the state updates should be outputted, 1 if the WiChR current should be outputted :return: states updates [1/ms] or WiChR current [pA/pF] """ # constants R = 8314.0 # mJ/mol/K T = 310.0 # K F = 96485.0 # C/mol if v != 0: FvRT = F * (v) / (R * T) else: FvRT = F * (1e-2) / (R * T) # ion concentrations nao = 110e-3 # M ko = 1e-3 # M # patch pipette parameters D_K = 187 # Diffusion coefficient in 10e-7 cm^2/s D_Na = 131 # Diffusion coefficient in 10e-7 cm^2/s tau_K = 6000 # ms tau_Na = D_K/D_Na * tau_K # ms na_pp = 1e-3 # M k_pp = 110e-3 # M # illumination opt_irr = 0 for opt_start in irr_times: if opt_start <= t < opt_start + opt_dur: opt_irr = opt_irr_max # states nai = X[0] # M ki = X[1] # M WiChR_C1 = X[2] WiChR_C1p = X[3] WiChR_O1 = X[4] WiChR_O2 = X[5] WiChR_O3 = X[6] WiChR_C2 = X[7] WiChR_O4 = X[8] if model == 'one-branch': # WiChR parameters PKPNa1 = 62.5 PKPNa2 = 52.4 WiChR_g2g1 = 0.904 WiChR_k = 0.229 # mm^2/mW/ms WiChR_a = 5.93 # 1/ms WiChR_e = 0.0033 # 1/ms/mV WiChR_d = 0.0035 # 1/ms # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R*T/F * np.log((PKPNa1*ko + nao)/(PKPNa1*ki + nai)) EK2 = R*T/F * np.log((PKPNa2*ko + nao)/(PKPNa2*ki + nai)) # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 WiChR_dC1pO1 = WiChR_a * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v+150) * WiChR_O1 WiChR_dO2C1 = WiChR_d * WiChR_O2 WiChR_dC1p = WiChR_dC1C1p - WiChR_dC1pO1 WiChR_dC1 = WiChR_dO2C1 - WiChR_dC1C1p WiChR_dO1 = WiChR_dC1pO1 - WiChR_dO1O2 WiChR_dO2 = WiChR_dO1O2 - WiChR_dO2C1 WiChR_dO3 = 0 WiChR_dC2 = 0 WiChR_dO4 = 0 # WiChR current I_WiChR1 = (WiChR_g * WiChR_O1) * (v - EK1) # pA/pF I_WiChR2 = (WiChR_g2g1 * WiChR_g * WiChR_O2) * (v - EK2) # pA/pF I_WiChR = I_WiChR1 + I_WiChR2 # pA/pF # Calculation of the K and Na current through WiChR (optional to use) WiChR_PNa1 = (I_WiChR1) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa1 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa2 = (I_WiChR2) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa2 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_I_K = FvRT * F * PKPNa1 * (WiChR_PNa1) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa2 * (WiChR_PNa2) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2) * (nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA elif model == 'two-branch': # WiChR parameters PKPNa1 = 79.1 PKPNa2 = 58.5 PKPNa3 = 54.1 PKPNa4 = 27.7 WiChR_g2g1 = 0.93 WiChR_g3g1 = 1.06 WiChR_g4g1 = 0.89 WiChR_k = 0.33 # mm^2/mW/ms if opt_irr == 0: WiChR_k1 = 0 else: WiChR_k1 = 24.4/(1+(50/opt_irr)**1.2) # 1/ms WiChR_k2 = 0.01 # mm^2/mW/ms WiChR_a = 6.0 # 1/ms WiChR_e = 0.0031 # 1/ms/mV WiChR_e1 = 0.0042 # 1/ms WiChR_d = 0.00315 # 1/ms WiChR_d1 = 0.0001 # 1/ms WiChR_d2 = 0.0044 # 1/ms WiChR_r = 0.00002 # 1/ms # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R * T / F * np.log((PKPNa1 * ko + nao) / (PKPNa1 * ki + nai)) EK2 = R * T / F * np.log((PKPNa2 * ko + nao) / (PKPNa2 * ki + nai)) EK3 = R * T / F * np.log((PKPNa3 * ko + nao) / (PKPNa3 * ki + nai)) EK4 = R * T / F * np.log((PKPNa4 * ko + nao) / (PKPNa4 * ki + nai)) # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 WiChR_dC1pO1 = WiChR_a * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v + 150) * WiChR_O1 WiChR_dO2C1 = WiChR_d * WiChR_O2 WiChR_dO2O3 = WiChR_k1 * WiChR_O2 WiChR_dO3O2 = WiChR_e1 * WiChR_O3 WiChR_dO3C2 = WiChR_d1 * WiChR_O3 WiChR_dC2O4 = WiChR_k2 * opt_irr * WiChR_C2 WiChR_dO4C2 = WiChR_d2 * WiChR_O4 WiChR_dC2C1 = WiChR_r * WiChR_C2 WiChR_dC1 = WiChR_dO2C1 + WiChR_dC2C1 - WiChR_dC1C1p WiChR_dC1p = WiChR_dC1C1p - WiChR_dC1pO1 WiChR_dO1 = WiChR_dC1pO1 - WiChR_dO1O2 WiChR_dO2 = WiChR_dO1O2 + WiChR_dO3O2 - WiChR_dO2C1 - WiChR_dO2O3 WiChR_dO3 = WiChR_dO2O3 - WiChR_dO3O2 - WiChR_dO3C2 WiChR_dC2 = WiChR_dO3C2 + WiChR_dO4C2 - WiChR_dC2O4 - WiChR_dC2C1 WiChR_dO4 = WiChR_dC2O4 - WiChR_dO4C2 # WiChR current I_WiChR1 = (WiChR_g * WiChR_O1) * (v - EK1) # pA/pF I_WiChR2 = (WiChR_g2g1 * WiChR_g * WiChR_O2) * (v - EK2) # pA/pF I_WiChR3 = (WiChR_g3g1 * WiChR_g * WiChR_O3) * (v - EK3) # pA/pF I_WiChR4 = (WiChR_g4g1 * WiChR_g * WiChR_O4) * (v - EK4) # pA/pF I_WiChR = I_WiChR1 + I_WiChR2 + I_WiChR3 + I_WiChR4 # pA/pF # Calculation of the K and Na current through WiChR (optional to use) WiChR_PNa1 = (I_WiChR1) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa1 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa2 = (I_WiChR2) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa2 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa3 = (I_WiChR3) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa3 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_PNa4 = (I_WiChR4) * 1 / (FvRT * F) * ((nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + PKPNa4 * (ki - ko * np.exp(-FvRT)) / ( 1 - np.exp(-FvRT))) ** (-1) WiChR_I_K = FvRT * F * PKPNa1 * (WiChR_PNa1) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa2 * (WiChR_PNa2) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa3 * (WiChR_PNa3) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) + \ FvRT * F * PKPNa4 * (WiChR_PNa4) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2 + WiChR_PNa3 + WiChR_PNa4) * (nai - nao * np.exp(-FvRT)) \ / (1 - np.exp(-FvRT)) # pA Q_K = -WiChR_I_K * C_cell # Total charge in pC n_current_K = Q_K / F # Moles of K due to current delta_concentration_current_K = n_current_K / V_cell / 1000 # Concentration change in M/ms Q_Na = -WiChR_I_Na * C_cell # Total charge in pC n_current_Na = Q_Na / F # Moles of Na due to current delta_concentration_current_Na = n_current_Na / V_cell / 1000 # Concentration change in M/ms # Update concentration due to diffusion delta_pipette_K = - 1 / tau_K * (ki - k_pp) # M/ms delta_pipette_Na = - 1 / tau_Na * (nai - na_pp) # M/ms dK = delta_concentration_current_K + delta_pipette_K # M/ms dNa = delta_concentration_current_Na + delta_pipette_Na # M/ms if flag_ode == 0: output = np.array([dNa, dK, WiChR_dC1, WiChR_dC1p, WiChR_dO1, WiChR_dO2, WiChR_dO3, WiChR_dC2, WiChR_dO4]) else: output = I_WiChR return output def simulate_WiChR(v, opt_irr_max, opt_dur, WiChR_g, V_cell, C_cell, model='one-branch', irr_times=[0], dt_sim=1, off_dur=2000, X0=[1e-3,110e-3,1,0,0,0,0,0,0]): """ Simulate the behavior of WiChR. :param v: transmembrane voltage in mV :param opt_irr_max: irradiance in mW/mm^2 :param opt_dur: duration of the light pulse in ms :param WiChR_g: WiChR conductance in nS/pF :param V_cell: cell volume in l :param C_cell: cell capacitance in F :param model: 'one-branch' or 'two-branch' photocycle model :param irr_times: list containing the illumination start times :param dt_sim: time step of the simulation in ms :param off_dur: duration after the light pulse that should be simulated :param X0: initial conditions of the WiChR model states :return: WiChR current in pA/pF """ time = np.arange(0, irr_times[-1] + opt_dur + off_dur, dt_sim) X = integrate.odeint(WiChR_model, X0, time, tfirst=True, args=(v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model), hmax=1) I_WiChR = [] for i in range(time.size): I_WiChR.append( WiChR_model(time[i], X[i, :], v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, flag_ode=1)) return time, np.array(I_WiChR)
