"""Contains opsin models and default parameters"""
from __future__ import annotations
from typing import Callable, Tuple
import warnings
from attrs import define, field, asdict, fields_dict
from brian2 import (
NeuronGroup,
Synapses,
Unit,
get_unit,
implementation,
check_units,
)
from nptyping import NDArray
from brian2.units import (
mm,
mm2,
nmeter,
Quantity,
second,
ms,
psiemens,
nsiemens,
mV,
volt,
amp,
mM,
)
from brian2.units.allunits import radian
from cleo.base import SynapseDevice
from cleo.light import LightDependent
# hacky MRO stuff: base class order is important here
[docs]@define(eq=False)
class Opsin(LightDependent, SynapseDevice):
"""Base class for opsin model.
Requires that the neuron model has a current term
(by default Iopto) which is assumed to be positive (unlike the
convention in many opsin modeling papers, where the current is
described as negative).
We approximate dynamics under multiple wavelengths using a weighted sum
of photon fluxes, where the ε factor indicates the activation
relative to the peak-sensitivy wavelength for an equivalent number of photons
(see Mager et al, 2018). This weighted sum is an approximation of a nonlinear
peak-non-peak wavelength relation; see ``notebooks/multi_wavelength_model.ipynb``
for details."""
@property
def action_spectrum(self):
"""Alias for ``light_receptor.spectrum``"""
return self.spectrum
@define(eq=False)
class MarkovOpsin(Opsin):
"""Base class for Markov state models à la Evans et al., 2016"""
required_vars: list[Tuple[str, Unit]] = field(
factory=lambda: [("Iopto", amp), ("v", volt)],
init=False,
)
@implementation(
"cython",
"""
cdef double f_unless_x0(double f, double x, double f_when_x0):
if x == 0:
return f_when_x0
else:
return f
""",
)
@check_units(f=1, x=volt, f_when_0=1, result=1)
def f_unless_x0(f, x, f_when_x0):
f[x == 0] = f_when_x0
return f
[docs]@define(eq=False)
class FourStateOpsin(MarkovOpsin):
"""4-state model from PyRhO (Evans et al. 2016).
rho_rel is channel density relative to standard model fit;
modifying it post-injection allows for heterogeneous opsin expression.
IOPTO_VAR_NAME and V_VAR_NAME are substituted on injection.
Defaults are for ChR2.
"""
g0: Quantity = 114000 * psiemens
gamma: Quantity = 0.00742
phim: Quantity = 2.33e17 / mm2 / second # *photon, not in Brian2
k1: Quantity = 4.15 / ms
k2: Quantity = 0.868 / ms
p: Quantity = 0.833
Gf0: Quantity = 0.0373 / ms
kf: Quantity = 0.0581 / ms
Gb0: Quantity = 0.0161 / ms
kb: Quantity = 0.063 / ms
q: Quantity = 1.94
Gd1: Quantity = 0.105 / ms
Gd2: Quantity = 0.0138 / ms
Gr0: Quantity = 0.00033 / ms
E: Quantity = 0 * mV
v0: Quantity = 43 * mV
v1: Quantity = 17.1 * mV
model: str = field(
init=False,
default="""
dC1/dt = Gd1*O1 + Gr0*C2 - Ga1*C1 : 1 (clock-driven)
dO1/dt = Ga1*C1 + Gb*O2 - (Gd1+Gf)*O1 : 1 (clock-driven)
dO2/dt = Ga2*C2 + Gf*O1 - (Gd2+Gb)*O2 : 1 (clock-driven)
C2 = 1 - C1 - O1 - O2 : 1
Theta = int(phi_pre > 0*phi_pre) : 1
Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1
Ga1 = k1*Hp : hertz
Ga2 = k2*Hp : hertz
Hq = Theta * phi_pre**q/(phi_pre**q + phim**q) : 1
Gf = kf*Hq + Gf0 : hertz
Gb = kb*Hq + Gb0 : hertz
fphi = O1 + gamma*O2 : 1
# v1/v0 when v-E == 0 via l'Hopital's rule
fv = f_unless_x0(
(1 - exp(-(V_VAR_NAME_post-E)/v0)) / ((V_VAR_NAME_post-E)/v1),
V_VAR_NAME_post - E,
v1/v0
) : 1
IOPTO_VAR_NAME_post = -g0*fphi*fv*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)
rho_rel : 1""",
)
extra_namespace: dict[str, Any] = field(
init=False, factory=lambda: {"f_unless_x0": f_unless_x0}
)
[docs] def init_syn_vars(self, opto_syn: Synapses) -> None:
for varname, value in {"C1": 1, "O1": 0, "O2": 0}.items():
setattr(opto_syn, varname, value)
[docs]@define(eq=False)
class BansalFourStateOpsin(MarkovOpsin):
"""4-state model from Bansal et al. 2020.
The difference from the PyRhO model is that there is no voltage dependence.
rho_rel is channel density relative to standard model fit;
modifying it post-injection allows for heterogeneous opsin expression.
IOPTO_VAR_NAME and V_VAR_NAME are substituted on injection.
"""
Gd1: Quantity = 0.066 / ms
Gd2: Quantity = 0.01 / ms
Gr0: Quantity = 3.33e-4 / ms
g0: Quantity = 3.2 * nsiemens
phim: Quantity = 1e16 / mm2 / second # *photon, not in Brian2
k1: Quantity = 0.4 / ms
k2: Quantity = 0.12 / ms
Gf0: Quantity = 0.018 / ms
Gb0: Quantity = 0.008 / ms
kf: Quantity = 0.01 / ms
kb: Quantity = 0.008 / ms
gamma: Quantity = 0.05
p: Quantity = 1
q: Quantity = 1
E: Quantity = 0 * mV
model: str = field(
init=False,
default="""
dC1/dt = Gd1*O1 + Gr0*C2 - Ga1*C1 : 1 (clock-driven)
dO1/dt = Ga1*C1 + Gb*O2 - (Gd1+Gf)*O1 : 1 (clock-driven)
dO2/dt = Ga2*C2 + Gf*O1 - (Gd2+Gb)*O2 : 1 (clock-driven)
C2 = 1 - C1 - O1 - O2 : 1
Theta = int(phi_pre > 0*phi_pre) : 1
Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1
Ga1 = k1*Hp : hertz
Ga2 = k2*Hp : hertz
Hq = Theta * phi_pre**q/(phi_pre**q + phim**q) : 1
Gf = kf*Hq + Gf0 : hertz
Gb = kb*Hq + Gb0 : hertz
fphi = O1 + gamma*O2 : 1
IOPTO_VAR_NAME_post = -g0*fphi*(V_VAR_NAME_post-E)*rho_rel : ampere (summed)
rho_rel : 1""",
)
[docs] def init_syn_vars(self, opto_syn: Synapses) -> None:
for varname, value in {"C1": 1, "O1": 0, "O2": 0}.items():
setattr(opto_syn, varname, value)
# TODO: not quite done
# @define(eq=False)
# class BansalThreeStatePump(MarkovOpsin):
# """3-state model from `Bansal et al. 2020 <10.1016/j.neuroscience.2020.09.022>`_.
# rho_rel is channel density relative to standard model fit;
# modifying it post-injection allows for heterogeneous opsin expression.
# IOPTO_VAR_NAME and V_VAR_NAME are substituted on injection.
# """
# Gd: Quantity = 0
# Gr: Quantity = 0
# ka: Quantity = 0
# p: Quantity = 0
# q: Quantity = 0
# phim: Quantity = 0
# E: Quantity = 0
# a: Quantity = 0
# b: Quantity = 0
# vartheta_max = 5 * mM
# kd = 16 * mM
# model: str = field(
# init=False,
# default="""
# dP0/dt = Gr*P6 - Ga*P0 : 1 (clock-driven)
# dP4/dt = Ga*P0 - Gd*P4 : 1 (clock-driven)
# P6 = 1 - P0 - P4 : 1
# Theta = int(phi_pre > 0*phi_pre) : 1
# Hp = Theta * phi_pre**p/(phi_pre**p + phim**p) : 1
# Ga = ka*Hp : hertz
# fphi = P4 : 1
# dCl_in/dt = a*(I_i + b*I_Cl_leak) : molar
# Cl_out : molar
# E_Cl = -26.67*mV * log(Cl_out/Cl_in) : volt
# I_Cl_leak = g_Cl * (E_Cl0 - E_Cl)
# Psi = vartheta_max*Cl_out / (kd + Cl_out) / 4.43 : 1
# I_i = fphi*(V_VAR_NAME_post-E)*Psi*rho_rel
# IOPTO_VAR_NAME_post = -(I_i + I_Cl_leak) : ampere (summed)
# rho_rel : 1""",
# )
# extra_namespace: dict[str, Any] = field(
# init=False, factory=lambda: {"E_Cl0": -70 * mV, "g_Cl": 2.3 * msiemens / cm2}
# )
# def init_opto_syn_vars(self, opto_syn: Synapses) -> None:
# raise NotImplementedError("Still need to figure out [Cl-_out]")
# opto_syn.P0 = 1
# opto_syn.P4 = 0
# opto_syn.P6 = 0
# opto_syn.Cl_out = 124 * mM
# opto_syn.Cl_in = np.exp(np.log(124) - 70 / 26.67) * mM
[docs]@define(eq=False)
class ProportionalCurrentOpsin(Opsin):
"""A simple model delivering current proportional to light intensity"""
I_per_Irr: Quantity = field(kw_only=True)
""" How much current (in amps or unitless, depending on neuron model)
to deliver per mW/mm2.
"""
# would be IOPTO_UNIT but that throws off Equation parsing
model: str = field(
init=False,
default="""
IOPTO_VAR_NAME_post = I_per_Irr / (mwatt / mm2)
* Irr_pre * rho_rel : IOPTO_UNIT (summed)
rho_rel : 1
""",
)
required_vars: list[Tuple[str, Unit]] = field(factory=list, init=False)
def __attrs_post_init__(self):
if isinstance(self.I_per_Irr, Quantity):
Iopto_unit = get_unit(self.I_per_Irr.dim)
else:
Iopto_unit = radian
self.per_ng_unit_replacements = [("IOPTO_UNIT", Iopto_unit.name)]
self.required_vars = [("Iopto", Iopto_unit)]