Merge branch 'greiner-reviewer-feedback' into 'main' (bb9c4ae1) · Commits · Public / WiChR_model

Readme.md

+2 −2

Original line number	Diff line number	Diff line
		@@ -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)`
		>

WiChR_model.py

+47 −44

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):
		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 = RT/F np.log((PKPNa1ko + nao)/(PKPNa1ki + 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)