change units of compute_current

jaxley_mech/channels/aoyama00.py

@@ -25,6 +25,7 @@ class Leak(Channel):
     """Leakage current"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -46,8 +47,8 @@ def compute_current(
     ):
         """Given channel states and voltage, return the current through the channel."""
         prefix = self._name
-        gLeak = params[f"{prefix}_gLeak"] * 1000  # mS/cm^2
-        return gLeak * (v - params[f"{prefix}_eLeak"])  # mS/cm^2 * mV = uA/cm^2
+        gLeak = params[f"{prefix}_gLeak"]  # S/cm^2
+        return gLeak * (v - params[f"{prefix}_eLeak"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -58,6 +59,7 @@ class Na(Channel):
     """Sodium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -89,11 +91,9 @@ def compute_current(
         """Given channel states and voltage, return the current through the channel."""
         prefix = self._name
         m, h = states[f"{prefix}_m"], states[f"{prefix}_h"]
-        gNa = (
-            params[f"{prefix}_gNa"] * (m**3) * h * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
+        gNa = params[f"{prefix}_gNa"] * (m**3) * h  # S/cm^2
 
-        return gNa * (v - params[f"{prefix}_eNa"])
+        return gNa * (v - params[f"{prefix}_eNa"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -124,6 +124,7 @@ class Kdr(Channel):
     """Delayed Rectifying Potassium Channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_params = {
             f"{self._name}_gKdr": 4.5e-3,  # S/cm^2
@@ -154,8 +155,8 @@ def compute_current(
         prefix = self._name
         m = states[f"{prefix}_m"]
         h = states[f"{prefix}_h"]
-        k_cond = params[f"{prefix}_gKdr"] * m**4 * h * 1000
-        return k_cond * (v - params["eK"])
+        k_cond = params[f"{prefix}_gKdr"] * m**4 * h  # S/cm^2
+        return k_cond * (v - params["eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -187,6 +188,7 @@ class Kto(Channel):
     """Transient Outward Potassium Channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_params = {
             f"{self._name}_gKto": 15e-3,  # S/cm^2
@@ -217,8 +219,8 @@ def compute_current(
         prefix = self._name
         m = states[f"{prefix}_m"]
         h = states[f"{prefix}_h"]
-        k_cond = params[f"{prefix}_gKto"] * m**3 * h * 1000
-        return k_cond * (v - params["eK"])
+        k_cond = params[f"{prefix}_gKto"] * m**3 * h  # S/cm^2
+        return k_cond * (v - params["eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -250,6 +252,7 @@ class Kar(Channel):
     """Anomalous rectifying Potassium Channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_params = {
             f"{self._name}_gKar": 4.5e-3,  # S/cm^2
@@ -276,8 +279,8 @@ def compute_current(
         """Compute the current through the channel."""
         prefix = self._name
         m = states[f"{prefix}_m"]
-        k_cond = params[f"{prefix}_gKar"] * m**5 * 1000
-        return k_cond * (v - params["eK"])
+        k_cond = params[f"{prefix}_gKar"] * m**5  # S/cm^2
+        return k_cond * (v - params["eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -299,6 +302,7 @@ class Ca(Channel):
     """L-type calcium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_params = {
             f"{self._name}_gCa": 9e-3,  # S/cm^2
@@ -330,8 +334,8 @@ def compute_current(
         """Compute the current through the channel."""
         prefix = self._name
         m = states[f"{prefix}_m"]
-        ca_cond = params[f"{prefix}_gCa"] * m**4 * 1000
-        current = ca_cond * (v - params["eCa"])
+        ca_cond = params[f"{prefix}_gCa"] * m**4  # S/cm^2
+        current = ca_cond * (v - params["eCa"])  # S/cm^2 * mV = mA/cm^2
         return current
 
     def init_state(self, states, v, params, delta_t):

jaxley_mech/channels/benison01.py

@@ -29,6 +29,7 @@ class Leak(Channel):
     """Leakage current"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -50,10 +51,8 @@ def compute_current(
     ):
         """""Return the updated states.""" ""
         prefix = self._name
-        gLeak = (
-            params[f"{prefix}_gLeak"] * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
-        return gLeak * (v - params[f"{prefix}_eLeak"])
+        gLeak = params[f"{prefix}_gLeak"]  # S/cm^2
+        return gLeak * (v - params[f"{prefix}_eLeak"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -64,6 +63,7 @@ class Na(Channel):
     """Sodium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -94,10 +94,8 @@ def compute_current(
         """Return the updated states."""
         prefix = self._name
         m, h = states[f"{prefix}_m"], states[f"{prefix}_h"]
-        gNa = (
-            params[f"{prefix}_gNa"] * (m**3) * h * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
-        return gNa * (v - params[f"{prefix}_eNa"])
+        gNa = params[f"{prefix}_gNa"] * (m**3) * h  # S/cm^2
+        return gNa * (v - params[f"{prefix}_eNa"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -128,6 +126,7 @@ class Kdr(Channel):
     """Delayed Rectifier Potassium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -145,7 +144,6 @@ def update_states(
         prefix = self._name
         m = states[f"{prefix}_m"]
         new_m = solve_gate_exponential(m, dt, *Kdr.m_gate(v))
-        # jax.debug.print("new_m={new_m}, v={v}", new_m=new_m, v=v)
         return {f"{prefix}_m": new_m}
 
     def compute_current(
@@ -154,11 +152,9 @@ def compute_current(
         """""Return the updated states.""" ""
         prefix = self._name
         m = states[f"{prefix}_m"]
-        gKdr = (
-            params[f"{prefix}_gKdr"] * (m**4) * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
+        gKdr = params[f"{prefix}_gKdr"] * m**4  # S/cm^2
 
-        return gKdr * (v - params[f"eK"])
+        return gKdr * (v - params[f"eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -186,6 +182,7 @@ class KA(Channel):
     """A type Potassium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -216,10 +213,8 @@ def compute_current(
         """Return the updated states."""
         prefix = self._name
         m, h = states[f"{prefix}_m"], states[f"{prefix}_h"]
-        gKA = (
-            params[f"{prefix}_gKA"] * (m**3) * h * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
-        return gKA * (v - params[f"eK"])
+        gKA = params[f"{prefix}_gKA"] * (m**3) * h  # S/cm^2
+        return gKA * (v - params[f"eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -250,6 +245,7 @@ class CaL(Channel):
     """L-type Calcium channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -280,10 +276,8 @@ def compute_current(
         """Return the updated states."""
         prefix = self._name
         m = states[f"{prefix}_m"]
-        gCaL = (
-            params[f"{prefix}_gCaL"] * (m**2) * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
-        current = gCaL * (v - states["eCa"])
+        gCaL = params[f"{prefix}_gCaL"] * (m**2)  # S/cm^2
+        current = gCaL * (v - states["eCa"])  # S/cm^2 * mV = mA/cm^2
         return current
 
     def init_state(self, states, v, params, delta_t):
@@ -306,6 +300,7 @@ class CaN(Channel):
     """N-type Ca channel"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -335,10 +330,8 @@ def compute_current(
         """Return the updated states."""
         prefix = self._name
         m, h = states[f"{prefix}_m"], states[f"{prefix}_h"]
-        gCaN = (
-            params[f"{prefix}_gCaN"] * (m**2) * h * 1000
-        )  # mS/cm^2, multiply with 1000 to convert Siemens to milli Siemens.
-        return gCaN * (v - states[f"eCa"])
+        gCaN = params[f"{prefix}_gCaN"] * (m**2) * h  # S/cm^2
+        return gCaN * (v - states[f"eCa"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -369,6 +362,7 @@ class CaPumpNS(Channel):
     """Non-spatial Calcium pump"""
 
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_constants = {
             "F": 96485.3329,  # C/mol (Faraday's constant)
@@ -406,7 +400,7 @@ def update_states(self, states, dt, v, params):
         K_pump = params["K_pump"]
         j_pump = v_pump * (Cai**4 / (Cai**4 + K_pump**4))
 
-        iCa = states["iCa"] / 1_000
+        iCa = states["iCa"]
 
         driving_channel = -10_000.0 * iCa / (2 * F * V_cell)
         driving_channel = select(
@@ -431,6 +425,7 @@ def init_state(self, states, v, params, delta_t):
 
 class KCa(Channel):
     def __init__(self, name: Optional[str] = None):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         prefix = self._name
         self.channel_params = {
@@ -464,9 +459,9 @@ def compute_current(
 
         # changed to the second-power from fm97
         x = (Cai / K_KCa) ** 2
-        gKCa = params[f"{prefix}_gKCa"] * (x / (1 + x)) * 1000
+        gKCa = params[f"{prefix}_gKCa"] * (x / (1 + x))  # S/cm^2
 
-        return gKCa * (v - params[f"eK"])
+        return gKCa * (v - params[f"eK"])  # S/cm^2 * mV = mA/cm^2
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state such at fixed point of gate dynamics."""
@@ -481,6 +476,7 @@ def __init__(
         self,
         name: Optional[str] = None,
     ):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         self.channel_constants = {
             "F": 96485.3329,  # C/mol (Faraday's constant)

jaxley_mech/channels/chen24.py

@@ -15,6 +15,7 @@
     "papers": [
         "Chen, Q., Ingram, N. T., Baudin, J., Angueyra, J. M., Sinha, R., & Rieke, F. (2024). Light-adaptation clamp: A tool to predictably manipulate photoreceptor light responses. https://doi.org/10.7554/eLife.93795.1",
     ],
+    "code": "https://github.com/chrischen2/photoreceptorLinearization",
 }
 
 
@@ -29,6 +30,7 @@ def __init__(
         atol: float = 1e-8,
         max_steps: int = 10,
     ):
+        self.current_is_in_mA_per_cm2 = True
         super().__init__(name)
         SolverExtension.__init__(self, solver, rtol, atol, max_steps)
         prefix = self._name
@@ -151,8 +153,13 @@ def compute_current(
             params[f"{prefix}_n"],
             params[f"{prefix}_k"],
         )
-        I = -k * G**n  # eq(4)
-        return I
+        I = -k * G**n  # eq(4) #pA
+
+        I *= 1e-9
+        area = 2 * jnp.pi * params["length"] * params["radius"] * 1e-8  # um^2 to cm^2
+        current_density = I / area
+
+        return current_density
 
     def init_state(self, states, v, params, delta_t):
         """Initialize the state at fixed point of gate dynamics."""