From 0899e4eb76c8a9f1fad83212ae9a5fc3084ad10e Mon Sep 17 00:00:00 2001 From: Joachim Greiner Date: Thu, 22 Jan 2026 12:26:34 +0000 Subject: [PATCH] Update model: Reviewer Feedback --- Readme.md | 4 +-- WiChR_model.py | 91 ++++++++++++++++++++++++++------------------------ 2 files changed, 49 insertions(+), 46 deletions(-) diff --git a/Readme.md b/Readme.md index faddac2..52062d4 100644 --- a/Readme.md +++ b/Readme.md @@ -18,7 +18,7 @@ Here we provide source code for the first experimentally validated photocycle mo To plot the WiChR current at 40 mV for 500 ms illumination with 4 mW/mm^2, given a whole-cell conductance of 4 nS/pF and a cell with 1.6 pl volume and 25.6 pF membrane capacitance, run: ->`time, I_WiChR = simulate_WiChR(40, 4, 500, 4, 1.6e-12, 25.6e-12)` +>`time, I_WiChR = simulate_WiChR(40, 4, 500, 4, 1.6e-12, 23.3e-12)` > >`plt.plot(time, I_WiChR)` > @@ -30,7 +30,7 @@ To plot the WiChR current at 40 mV for 500 ms illumination with 4 mW/mm^2, given If the photocurrent should be calculated with the two-branch model, use the option model='two-branch': ->`time, I_WiChR = simulate_WiChR(40, 4, 500, 4, 1.6e-12, 25.6e-12, model='two-branch')` +>`time, I_WiChR = simulate_WiChR(40, 4, 500, 4, 1.6e-12, 23.3e-12, model='two-branch')` > >`plt.plot(time, I_WiChR)` > diff --git a/WiChR_model.py b/WiChR_model.py index fa1a574..2062d67 100644 --- a/WiChR_model.py +++ b/WiChR_model.py @@ -1,10 +1,9 @@ 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): +def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, T, 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 @@ -16,13 +15,14 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma :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 T: Temperature in Kelvin :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 + deltaT = T - 296 F = 96485.0 # C/mol if v != 0: FvRT = F * (v) / (R * T) @@ -34,8 +34,8 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma 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 + D_K = 196 # Diffusion coefficient in 10e-7 cm^2/s + D_Na = 133 # 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 @@ -60,13 +60,15 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma 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 + PKPNa1 = 60.0 + PKPNa2 = 52.5 + WiChR_g2g1 = 0.91 + WiChR_k = 0.238 # mm^2/mW/ms + WiChR_a = 0.0351 # 1/ms/mV + WiChR_e = 0.0034 # 1/ms/mV + WiChR_d = 0.0034 # 1/ms + WiChR_Q10d = 1.52 + # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R*T/F * np.log((PKPNa1*ko + nao)/(PKPNa1*ki + nai)) @@ -74,9 +76,9 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 - WiChR_dC1pO1 = WiChR_a * WiChR_C1p + WiChR_dC1pO1 = WiChR_a * (150-v) * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v+150) * WiChR_O1 - WiChR_dO2C1 = WiChR_d * WiChR_O2 + WiChR_dO2C1 = WiChR_d * WiChR_Q10d**(deltaT/10) * WiChR_O2 WiChR_dC1p = WiChR_dC1C1p - WiChR_dC1pO1 WiChR_dC1 = WiChR_dO2C1 - WiChR_dC1C1p @@ -101,31 +103,29 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma 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 + FvRT * F * PKPNa2 * (WiChR_PNa2) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA/pF + WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2) * (nai - nao * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA/pF 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 + PKPNa1 = 64 + PKPNa2 = 52 + PKPNa3 = 43 + PKPNa4 = 28 + WiChR_g2g1 = 0.91 + WiChR_g3g1 = 1.1 + WiChR_g4g1 = 1.0 + WiChR_k = 0.28 # mm^2/mW/ms + WiChR_k1 = 0.097 # mm^2/mW/ms + WiChR_k2 = 0.016 # mm^2/mW/ms + WiChR_a = 0.035 # 1/ms + WiChR_e = 0.0032 # 1/ms/mV + WiChR_e1 = 0.004 # 1/ms + WiChR_d = 0.0034 # 1/ms + WiChR_d1 = 0.00022 # 1/ms + WiChR_d2 = 0.004 # 1/ms + WiChR_r = 0.00022 # 1/ms + WiChR_Q10d = 1.5 # Calculation of the reversal potential based on the selectivity ratio PK/PNa EK1 = R * T / F * np.log((PKPNa1 * ko + nao) / (PKPNa1 * ki + nai)) @@ -135,10 +135,10 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma # Updates of the states WiChR_dC1C1p = WiChR_k * opt_irr * WiChR_C1 - WiChR_dC1pO1 = WiChR_a * WiChR_C1p + WiChR_dC1pO1 = WiChR_a * (150 - v) * WiChR_C1p WiChR_dO1O2 = WiChR_e * (v + 150) * WiChR_O1 - WiChR_dO2C1 = WiChR_d * WiChR_O2 - WiChR_dO2O3 = WiChR_k1 * WiChR_O2 + WiChR_dO2C1 = WiChR_d * WiChR_Q10d**(deltaT/10) * WiChR_O2 + WiChR_dO2O3 = WiChR_k1 * opt_irr * WiChR_O2 WiChR_dO3O2 = WiChR_e1 * WiChR_O3 WiChR_dO3C2 = WiChR_d1 * WiChR_O3 WiChR_dC2O4 = WiChR_k2 * opt_irr * WiChR_C2 @@ -178,9 +178,9 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma 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 + FvRT * F * PKPNa4 * (WiChR_PNa4) * (ki - ko * np.exp(-FvRT)) / (1 - np.exp(-FvRT)) # pA/pF WiChR_I_Na = FvRT * F * (WiChR_PNa1 + WiChR_PNa2 + WiChR_PNa3 + WiChR_PNa4) * (nai - nao * np.exp(-FvRT)) \ - / (1 - np.exp(-FvRT)) # pA + / (1 - np.exp(-FvRT)) # pA/pF Q_K = -WiChR_I_K * C_cell # Total charge in pC n_current_K = Q_K / F # Moles of K due to current @@ -207,7 +207,7 @@ def WiChR_model(t, X, v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_ma -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]): +def simulate_WiChR(v, opt_irr_max, opt_dur, WiChR_g, V_cell, C_cell, model='one-branch', T=310, 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 @@ -217,6 +217,7 @@ def simulate_WiChR(v, opt_irr_max, opt_dur, WiChR_g, V_cell, C_cell, model='one- :param V_cell: cell volume in l :param C_cell: cell capacitance in F :param model: 'one-branch' or 'two-branch' photocycle model + :param T: Temperature in Kelvin :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 @@ -225,12 +226,14 @@ def simulate_WiChR(v, opt_irr_max, opt_dur, WiChR_g, V_cell, C_cell, model='one- """ 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) + 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, T), 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)) + WiChR_model(time[i], X[i, :], v, WiChR_g, V_cell, C_cell, irr_times, opt_dur, opt_irr_max, model, T, flag_ode=1)) return time, np.array(I_WiChR) + + -- GitLab