Double Two-State Opsin Model With Autonomous Parameter Inference

1. Introduction

With optogenetics, neuronal firing can be controlled with light. This is achieved by genetically expressing opsins, light sensitive ion channels or pumps, in cells or cell subtypes. The merger of this genetic expression and optical stimulation results in superior spatiotemporal resolution with respect to the conventional neuromodulation techniques. Consequently, it is an ideal investigative tool for behavioral studies and a promising biomedical treatment for medical disorders such as epilepsy, Parkinson's disease and beyond the brain conditions (Aravanis et al., 2007; Abilez et al., 2011; Gerits and Vanduffel, 2013; Williams et al., 2013; Klapoetke et al., 2014; Carrette et al., 2015; Chen et al., 2015; Tønnesen and Kokaia, 2017).

The first light sensitive ion channels were discovered in the green alga Chlamydomonas reinhardtii by Nagel et al. (2002). Genetic engineering has led to a variety of opsins, such as red-shifted, step-function and ultrafast opsins, and mutants with altered ion selectivity (Gunaydin et al., 2010; Gerits and Vanduffel, 2013; Azimihashemi et al., 2014). An example of the latter is ChR2(H134R), which is addressed in this paper. Furthermore, other natural versions are continuously being discovered as well. An example are the MerMAIDs, which is a family of metagenomically discovered marine anion-conducting and intensely desensitizing channelrhodopsins (Oppermann et al., 2019).

In its initial dark-adapted (IDA) state and under voltage clamp conditions, ChR2's photocurrent exhibits a peak (I_peak) and a steady-state current (I_ss) (Bruun et al., 2015). The peak is reached within 1-2 ms and followed by fast decay onto a steady-state plateau. This is due to light adaptation (Figure 1A, left). Post-illumination, there is a bi-exponential decay back to baseline, rendering the channel in an apparent dark-adapted state (DA_app). This is observed by applying a second stimulation after a short period of time (<10 s), which results in a reduced transient response with a maintained steady-state current (Figure 1A, right) (Nikolic et al., 2009; Bruun et al., 2015).

FIGURE 1

Figure 1. The Channelrhodopsin-2 photocurrent and photocycle. (A), The photocurrent for a single light pulse on the left. Right, response to a S1-S2 pulse protocol with variable inter-pulse intervals. Light pulses are indicated with blue bars and target features with black arrows. (B), The unified photocycle as proposed by Kuhne et al. (2019). (C), Previously proposed models. (C, top left) a three state cycle model with second light dependent step (dotted or dashed step) (Ernst et al., 2008). (C top right) a four state branching model (Williams et al., 2013). (C bottom) a six state model with two extra activation intermediates (Grossman et al., 2013). (D), The proposed double two-state opsin model (22OM) with separation of open-closing mechanism and conductance change due to dark-light adaptation. The latter is captured in the mathematical R and S model state pair. DA and LA indicate dark and light adapted molecule states, respectively. O means open, C is closed and D is desensitized. Blue arrows indicate light dependent rates.

ChR2 comprises seven transmembrane helices. These are covalently bound with a retinal chromophore forming a protonated retinal Schiff base (RSBH⁺). In its IDA (D470), retinal is in an all-trans configuration (Bruun et al., 2015). Upon illumination, a 13 trans-cis isomerization is triggered that initiates a cascade of conformational changes with opening of the pore as result (P520). Before returning back to the dark adapted state, the channel converts to a non-conducting state P470. This happens on a millisecond timescale, while complete recovery takes seconds (Stehfest and Hegemann, 2010; Ritter et al., 2013; Volkov et al., 2017). There is strong evidence for a second photocycle, with similar intermediates (Stehfest and Hegemann, 2010; Schneider et al., 2013; Bruun et al., 2015; Deisseroth and Hegemann, 2017). However, the transition between the two photocycles is still under debate (Bruun et al., 2015; Deisseroth and Hegemann, 2017). Recently, Kuhne et al. (2019) proposed an unifying ChR2 photocycle model consisting of two parallel photocycles, with three reaction pathways (Figure 1B).

In silico, the photocurrent is currently modeled with either a three- or a four-state Markov model (Figure 1C). This is in accordance with the single and double photocycle hypothesis, respectively. The opening is reduced to a single state transition. This is because the D480 → P500 and P500 → P390 transitions occur on a much faster timescale. However, in order to represent fast closure, slow recovery and a steady-state current, a second photon absorption step is proposed for the three state model (Nagel et al., 2003; Stehfest and Hegemann, 2010). The photochemical transition either increases the recovery rate or acts as equilibrium modulator between the open and desensitized state. The six state model, as depicted in Figure 1C, bottom, is an extended version of the four state model. The additional two intermediates are to correctly account for the activation time after retinal isomerizations and to avoid explicit time dependent rates (Grossman et al., 2013). The four state model is in agreement with the second photocycle hypothesis with modeling of two open and closed states. The transition as depicted in Figure 1C, top right is according to the older transition hypothesis, not to the latest unifying photocycle model of Kuhne et al. (2019).

In silico studies, which allow for extensive and systematical investigation of the effects of the current kinetics, require an accurate model of the to be investigated opsin. To date, an accurate model consists of four differential equations (Williams et al., 2013). Using such a model therefore increases the computational burden enormously, especially in case of multi-compartment or network studies. Moreover, due to the expanding possibilities, selection of the correct opsin is crucial to have an optimal optogenetic tool. These four state Markov models are not easily fit as they require preliminary knowledge of the parameter space and its complex interactions. Furthermore, finding the optimal parameters is time-consuming as the set of differential equations has to be evaluated at each step in the selected optimization algorithm.

In this study we propose, for the first time (to the authors' knowledge), the use of a double two-state model structure for modeling of the opsin's photocurrents (Figure 1D). A fit is created of the ChR2(H134R) mutant and compared to the 4SB model of Williams et al. (2013). The performance of both models are tested in a regular spiking neuron (Pospischil et al., 2008). The difference in computation speed is assessed as well, this in the aforementioned regular spiking neuron for different stimulation patterns and in the sparse Pyramidal-Interneuron-Network-Gamma (sPING) network model (Börgers and Kopell, 2005; Sherfey et al., 2018) with increasing number of transfected neurons. Finally, the versatility of the proposed model is evaluated with a fit to a MerMAID opsin. The improvements of this work with respect to the current state of the art are:

• A double two-state opsin model structure resulting in a reduced complexity of fifty percent, which leads to an increase of computation speed valuable for optogenetic neurostimulation in conductance based models.

• A fitting procedure that allows for autonomous and accurate model fitting starting from a vast parameter space. Moreover, the fit time is reduced significantly by using an analytical solution to the set of differential equations (describing the double two-state model structure), possible under voltage clamp conditions and rectangular optical pulses.

• The proposed model structure can be used to accurately describe distinctive opsins.

However, the model does not include the non-instantaneous response of the retinal complex to light. Consequently, the response to short (< τ_on) light pulses might be overestimated. Preliminary results of this work have been presented on the 28th Annual Computational Neuroscience Meeting (CNS*2019) (Schoeters et al., 2019).

2. Materials and Methods

In this study, a double two-state opsin model structure (22OM) was tested as alternative for opsin modeling. Below, we first describe the model in full and indicate the link between parameters and certain features. Next, the fitting procedure is elaborated. Finally, we describe the models and metrics used in the analysis of the model performance and computational speed.

2.1. The Model

The proposed model is based on the original voltage gated sodium model of Hodgkin and Huxley (Hodgkin and Huxley, 1952). It consists of two independent two-state pairs as depicted in Figure 1D. In contrast to the sodium model, where the second two-state pair represents the inactivation gate, it represents here the change in conductance due to dark-light adaptation.

After a long enough dark period, the molecules are assumed to be all in closed, dark adapted state. Upon stimulation, the channel opens with a transition C → O. On a slightly slower time scale the equilibrium between dark and light adapted molecules is reached. Light adapted molecules have a lower conductance than those that are in the dark adapted state. This change in conductance is captured by a transition R → S. The relationship between these mathematical model states and the physical dark and light adapted states of the opsin molecules is obtained via a linear transformation, i.e., R = (g_ChR2·DA + g_LA·LA)/g_ChR2. Consequently, R (S) is one (zero) when fully dark adapted and g_LA/g_ChR2 (respectively, 1 − g_LA/g_ChR2) when fully light adapted, with g_LA the conductivity of a light adapted channel. DA and LA are the possibilities of the opsin molecules being in a dark or light adapted state, respectively. By using the R state in the model, g_LA does not need to be determined, therefore reducing the number of model parameters. The established equilibria of both state pairs depend on the level of optical excitation. After photostimulation, the channels close (O → C). Moreover, they all return to the dark adapted state after a long enough recovery period, which is on a much slower time scale than the other temporal kinetics. Because of this slower time scale, the transition S → R has to be light dependent as well. Otherwise the equilibrium would be completely on the side of S for every optical excitation level. The ChR2 photocurrent can thus be determined as follows:

$\begin{array}{ll}{i}_{\text{ChR}2}=g_{\text{ChR}2}G(V)(O·R)(V-E_{\text{ChR}2})& (1)\end{array}$

with

$\begin{array}{ll}\frac{dO}{dt}=\frac{O_{\infty }(I,V)-O(t)}{{\tau }_O(I,V)}& (2)\end{array}$

$\begin{array}{ll}\frac{dR}{dt}=\frac{R_{\infty }(I,V)-R(t)}{{\tau }_R(I,V)}& (3)\end{array}$

where g_ChR2 is the maximal specific conductivity of the fully dark adapted channel, G(V) is a rectification function, V the membrane potential, I the light intensity, E_ChR2 the equilibrium potential and O the fraction of molecules in the open state, with O_∞ and τ_O its corresponding equilibrium and time constant. R_∞ and τ_R are the respective equilibrium and time constants of the R state.

Under voltage clamp conditions and a rectangular optical pulse with constant light intensity, the photocurrent can be expressed in a closed form analytical expression:

$\begin{array}{l}{i}_{\text{ChR2}}=g_{\text{ChR2}}G(V)(O_{\text{ChR2}}^{\text{on}}(t)+O_{\text{ChR2}}^{\text{off}}(t))·(R_{\text{ChR2}}^{\text{on}}(t)\\ +R_{\text{ChR2}}^{\text{off}}(t))(V-E_{\text{ChR2}})& (4)\end{array}$

with

$\begin{array}{l}{O}_{ChR2}^{\text{on}}(t)=\left[O_{\text{∞}}-(O_{\text{∞}}-O_0)exp\left(-\frac{t-t_{\text{on}}}{{\tau }_{\text{O}}(I,V)}\right)\right]\\ ·\Theta (t-t_{\text{on}})\Theta (t_{\text{off}}-t)& (5)\end{array}$

$\begin{array}{ll}{O}_{\text{ChR}2}^{\text{off}}(t)=O_{\text{ChR}2}^{\text{on}}(t_{\text{off}})exp\left(-\frac{t-t_{\text{off}}}{{\tau }_{\text{O}}(0,V)}\right)\Theta (t-t_{\text{off}})& (6)\end{array}$

$\begin{array}{l}{R}_{\text{ChR}2}^{\text{on}}(t)=\left[R_{\text{∞}}-(R_{\text{∞}}-R_0)exp\left(-\frac{t-t_{\text{on}}}{{\tau }_{\text{R}}(I,V)}\right)\right]\\ ·\Theta (t-t_{\text{on}})\Theta (t_{\text{off}}-t)& (7)\end{array}$

$\begin{array}{l}{R}_{\text{ChR}2}^{\text{off}}(t)=\left[1-\left(1-R_{\text{ChR}2}^{on}(t_{\text{off}})\right)·exp\left(-\frac{t-t_{\text{off}}}{{\tau }_{\text{R}}(0,V)}\right)\right]\\ \Theta (t-t_{\text{off}})& (8)\end{array}$

with Θ the Heaviside function, O₀ and R₀ the initial values of O and R at t = t_on (respectively, 0 and 1 when fully dark adapted) and, t_on and t_off, respectively, the onset and offset of the optical pulse.

This is of particular use during the fitting procedure as the model is fit to experimental data, recorded under the same aforementioned conditions. Moreover, strong correlations between the model time constants and experimentally determined features (Figure 1A) are observed. These can be exploited to obtain a first approximation of the model's parameters (see section 2.2). When τ_O ≪ τ_R, the transition rate time constant τ_O can be easily obtained from the activation (τ_on) and deactivation (τ_off) time constants. Under the same conditions, τ_R strongly correlates with the inactivation time constant (τ_inact) when I ≠ 0. The recovery time constant needs to be scaled as shown in Equation 10 to get a good approximation of the dark-light adaptation time constant under dark (I = 0) conditions. This relationship is obtained by evaluating the recovery time definition with the given model equations, i.e., τ_recov = t_on,2 − t_off,1 → I_p,2/I_p,1 = 1 − exp(−1). Here, t_on,2 is the onset time of the second pulse, t_off,1 the offset of the first pulse, and I_p,2 and I_p,1 the current peak value of second and first pulse, respectively.

$\begin{array}{ll}{\tau }_{\text{O}}(I,V)\approx {\tau }_{\text{on}},{\tau }_{\text{O}}(0,V)\approx {\tau }_{\text{off}}\text{ and }{\tau }_{\text{R}}(I,V)\approx {\tau }_{\text{inact}}& (9)\end{array}$

$\begin{array}{ll}{\tau }_{\text{R}}(0,V)\approx {\tau }_{\text{recov}}/\left(1-ln\frac{1}{1-I_{\text{ratio}}}\right)& (10)\end{array}$

Furthermore, following conditions need to be met for the relationship to hold true:

$\begin{array}{l}{t}_{\text{p},1}-t_{\text{on},1}\approx t_{\text{p},2}-t_{\text{on},2}\\ t_{\text{p},1}-t_{\text{on},1}>{\tau }_{\text{O}}\\ t_{\text{off},1}-t_{\text{on},1}>{\tau }_{\text{R}}& (11)\end{array}$

The first, t_p,i − t_on,i is the time required to reach the peak value since onset of pulse i. This needs to be approximately the same in both first and second pulse, while these need to be significantly larger than the activation time constant. The last one requires that the steady-state value is reached at the end of the first pulse.

Unless specified, the time constants and time in this study are in seconds, the membrane potential in mV and the intensity in W/m². The units of the conductance depend on the experimental data of each opsin, i.e., mS/cm² and μS in case of the ChR2(H134R) and MerMAID fit, respectively.

2.2. The Fitting Procedure

Due to the dependency on both the potential and light intensity, more than twenty parameters need to be inferred. This vast parameter space impedes finding the optimal solution which is at a high computational cost. To alleviate this, the fitting procedure can be divided into four steps.

The first step is the extraction of the features, which is described by Williams et al. (2013). The peak current (I_peak) is the maximal deflection from baseline. The steady-state current (I_ss) is the plateau value. The current ratio (I_ratio) is then I_ss/I_peak. The time constants are extracted using mono-exponential curve fits. To this end, a nonlinear least-squares curve fit is performed, with a trust-region-reflective algorithm. Furthermore, a multi-start algorithm with ten starting points was used to ensure finding of the global solution. The variable and function tolerance were set to 10⁻¹². The recovery time constant, i.e., the time necessary between two pulses to have a second peak current which is 63% of the first peak (see definition in previous subsection), was determined from a set of two-pulse experiments.

Next, τ_O and τ_R are fit to the obtained target data. Both are fit to the corresponding time constants (see Equations 9, 10) using the aforementioned nonlinear least-squares method. Again, a multi-start algorithm is used but with 2000 starting points. For the intensity dependence, sigmoidal functions on the log-scale are used while for the voltage dependence a logistic regression was selected. The two dependencies are combined by either a multiplication or a reciprocal addition. The relationships and combination schemes are given by Equations 12–16, with p_i, i = 1 → 6 indicating the unknown parameters of each relationship individually.

$\begin{array}{ll}{\tau }_{\text{O}}(I)=\frac{p_3}{1+exp(p_1/p_2)·I^{1/p_2·ln(10)}}& (12)\end{array}$

$\begin{array}{l}{\tau }_{\text{R}}(I)=p_1(1-\frac{p_2}{1+exp(p_3/p_4)·I^{-1/p_4·ln(10)}}\\ {\displaystyle -\frac{(1-p_2)}{1+exp(p_5/p_6)·I^{-1/p_6·ln(10)}})}& (13)\end{array}$

$\begin{array}{ll}{\tau }_{\text{X}}(V)=\frac{p_1}{1+exp(-(V-p_2)/p_3)}& (14)\end{array}$

$\begin{array}{lll}{\tau }_{\text{X}}(I,V)=& {\tau }_{\text{X}}(I)·{\tau }_{\text{X}}(V)& (15)\end{array}$

or

$\begin{array}{ll}{\left[{\left({\tau }_{\text{X}}(I)\right)}^{-1}+{\left({\tau }_{\text{X}}(V)\right)}^{-1}\right]}^{-1}& (16)\end{array}$

$\begin{array}{ll}{O}_{\text{∞}}(I)=\frac{1}{1+exp(p_1/p_2)·I^{-1/p_2·ln(10)}}& (17)\end{array}$

$\begin{array}{ll}{R}_{\text{∞}}(I)=1-\frac{p_3}{1+exp(p_1/p_2)·I^{-1/p_2·ln(10)}}& (18)\end{array}$

$\begin{array}{ll}G(V)=\frac{p_1·\left(1-p_{2}exp(-(V-E_{\text{ChR}2})/p_3)\right)}{V-E_{\text{ChR}2}}& (19)\end{array}$

In a third step, the parameters of the rectification function G(V) and the equilibrium constants O_∞ and R_∞ are fit. The used relationships are given in Equations 17–19, respectively. The potential dependence of O_∞ and R_∞ are omitted because this is mostly covered by the rectification function. The parameter values are determined by minimizing the cost function described below:

$\begin{array}{l}{f}_{\text{cost}}={\left(\frac{1}{N}\left[{\displaystyle \sum _{i=1\to N}}\Delta I_{\text{peak}}{(I_i,V_i)}^2+\Delta I_{\text{ss}}{(I_i,V_i)}^2+\Delta I_{\text{ratio}}{(I_i,V_i)}^2\right]\right)}^{1/2}\\ \Delta I_{\text{x}}(I_i,V_i)=w_{\text{x}}(y_{\text{x}}(I_i,V_i)-t_{\text{x},{\text{I}}_{\text{i}},{\text{V}}_{\text{i}}}),\text{with }x=peak,ss,ratio& (20)\end{array}$

Here, y_x and t_x,_Ii,Vi are, respectively the model output and target value at stimulation values (I,V), with $y_{\text{peak}}={max}_{\text{t}}(|i_{\text{ChR}2}^{\text{on}}(t,I,V)|)$, $y_{\text{ss}}=i_{\text{ChR}2}^{\text{on}}(t_{\text{off}})$ and y_ratio = y_ss/y_peak. $i_{\text{ChR}2}^{\text{on}}(t,I,V)$ is the current during the photostimulation pulse (t ϵ [t_on, t_off]) for a certain intensity I and voltage V. The current is calculated by evaluating Equations 4–8 with the determined dependencies in the previous step. N is the total number of stimulation sets (I,V). The minimization of f_cost is performed with the MATLAB fmincon-function and multi-start algorithm with 3000 starting points to increase chance of finding the global optimum. The upper and lower boundaries as well as the initial conditions are summarized in Table 1. Extra nonlinear constraints are applied to assure that O_∞ approaches one for high intensities (see sections 3.1, 3.4) and G(V) ≥ 0. A final constraint ensures a current decay back to baseline after the optical stimulation, i.e., i_on(t_off) > i_off(t) or $O^{\text{on}}(t_{\text{off}})·R^{\text{on}}(t_{\text{off}})>O^{\text{off}}(t)·R^{\text{off}}(t)$, resulting in:

$\begin{array}{ll}{R}_{\text{∞}}(I,V)>1-\frac{{\tau }_{\text{R}}(0,V)}{{\tau }_{\text{R}}(0,V)+{\tau }_{\text{O}}(0,V)}& (21)\end{array}$

Finally, a global optimization is performed with the parameters of all rate functions included. First, a new parameter space is defined, which is 10% of the original parameter space but centered around the values obtained in previous steps and limited by the former. With the gathered dependencies, the ChR2 current is calculated according to Equation 4. All model features are now extracted in the same manner as performed on the experimental data. These are used to determine a cost function which is the weighted root mean square error Equation 20, with additional terms: $\Delta {\tau }_{\text{on}}{(I,V)}^2$, $\Delta {\tau }_{\text{off}}{(I,V)}^2$, $\Delta {\tau }_{\text{inact}}{(I,V)}^2$ and $\Delta {\tau }_{\text{rec}}{(I,V)}^2$. Subsequently, the problem is optimized with a bounded particle swarm optimization (Hassan et al., 2005; Helwig and Wanka, 2007; Chen, 2018), containing 1000 particles and with a time limit of 24 h. The same solver settings and constraints are imposed as described in previous steps. The single-pulse experiments are evaluated with a time step of 1.5·10⁻⁴s, while for the two-pulse experiments a step of 1 ms is used.

TABLE 1

Table 1. Summary of all parameters.

2.3. Performance Tests

In this study, two opsin fits were performed. First, a fit is made to the data reported by Williams et al. (2013) of the ChR2(H134R) (Williams et al., 2013). The model accuracy is compared to the four state Markov model created by the same group. Four metrics are used to analyze the goodness-of-fit, i.e., Root mean square error (RMSE), Root mean square normalized error (RMSNE), Root mean square weighted error (RMSWE) and root mean square Z-score error (RMSZE):

$\begin{array}{ll}RMSWE={\left(\frac{1}{N}{\displaystyle \sum _{i=1\to N}}{w}_x^2·{\left[y_x(I_i,V_i)-t_{x,I_i,V_i}\right]}^2\right)}^{1/2}& (22)\end{array}$

where w_x equals 1, 1/t_x,_Ii,Vi, or 1/σ_x,_Ii,Vi in case of RMSE, RMSNE or RMSZE, respectively. y_x(I_i, V_i), t_x,_Ii,Vi and σ_x,_Ii,Vi are the model output, target feature and standard deviation of target feature x under intensity I and voltage V of set i, and w_x are the weights used in f_cost. The metrics are also determined in the overall, time constant features (τ_on + τ_off + τ_inact + τ_rec) only and current features (I_p + I_ss + I_ratio) only case. Here the squared errors of all features are summed first before taking the root and mean. The RMSWE is equivalent to the training error. However, it could not be used to compare the model fits as the used weights were not equal across fitting procedures (different weights were used in the 4SB fit, see Williams et al., 2013). Therefore, the other metrics were defined as well. Where the RMSE is biased by high values, the RMSNE is biased by values close to zero and RMSZSE which includes the uncertainty of the target features via σ_x,_Ii,Vi but could not be determined for the recovery time constant.

Both models are then implemented in a regular spiking neuron, described in Pospischil et al. (2008). The strength duration curves (SDC) are determined. When the irradiance is selected as strength for the SDC, a poor fit is obtained. This is due to the assumption of an RC equivalent circuit and a rectangular stimulation pulse in the Hill-Lapicque relationship Equation 23 (Noble and Stein, 1966; Williams and Entcheva, 2015). Therefore, the SDC fit is performed on the average inward stimulation current or temporal averaged current (i_ChR2,avg, TAC), as described by Williams and Entcheva (2015).

$\begin{array}{ll}{i}_{\text{ChR}2,\text{avg}}=\frac{I_{\text{ChR}2,\text{rheo}}}{\left(1-exp(-\frac{PD}{{\tau }_{\text{ChR}2,\text{chron}}/ln(2)})\right)}& (23)\end{array}$

$\begin{array}{ll}{i}_{\text{ChR}2,\text{avg}}=\frac{1}{PD}·{\displaystyle {\int }_0^{T_{\text{end}}}}{i}_{\text{ChR}2}(t)dt& (24)\end{array}$

with PD the pulse duration and T_end one second after the end of the pulse. The relationship between the irradiance and i_ChR2,avg is obtained through a power series fit, which allows calculation of the irradiance rheobase (I_rheo) and chronaxie (τ_chron) as follows:

$\begin{array}{ll}{I}_{\text{rheo}}=a·I_{\text{ChR}2,\text{rheo}}^b+c& (25)\end{array}$

$\begin{array}{ll}{\tau }_{\text{chron}}=-\frac{{\tau }_{\text{ChR}2,\text{chron}}}{ln(2)}ln\left(1-\frac{I_{\text{ChR}2,\text{rheo}}}{{\left[(2I_{\text{rheo}}-c)/a\right]}^{1/b}}\right)& (26)\end{array}$

where a, b and c are parameters obtained in an empirically power series fit of the irradiance curve vs. the inward stimulation current ($I=a·{(i_{\text{ChR}2,\text{avg}})}^b+c$) (Williams and Entcheva, 2015).

Moreover, the simulation speed is determined for different stimulation paradigms, i.e., simulation time (T_end)/runtime in a regular spiking neuron (Pospischil et al., 2008). Therefore, we varied the pulse repetition frequency, stimulation time and duty cycle. The intensity was fixed for each model and set to a value that elicited a firing rate of 100 Hz in the regular spiking neuron in case of a two pulse stimulation of 2 s with duty cycle 0.5 and pulse repetition frequency of 1 Hz. The models were solved by the MATLAB Variable Step Variable Order solver (VSVO) ode113-solver (order 1-13, Adams-Bashort-Moulton predictor-corrector pairs) (Shampine and Reichelt, 1997), with a maximum time step of 100 μs and default tolerances, i.e., relative and absolute tolerance equal to 10⁻³ and 10⁻⁶, respectively.

Finally, computational gain with the proposed model compared to the 4 state Markov model was tested in a network model with an increasing number of transfected neurons. Therefore, we used the sparse Pyramidal-Interneuron-Network-Gamma (sPING) (Börgers and Kopell, 2005), which was implemented via the DynaSim toolbox (Sherfey et al., 2018). The ChR2(H134R) models were added to the pyramidal neurons. The number of inhibitory neurons was varied between 3 and 100 while the 4/1, pyramidal/interneuron-ratio was maintained. The network was fully connected and the GABAa and AMPA conductivities were scaled such that the total input per neuron stayed the same, i.e., $g_{\text{GABAa}}=2/(N_{\text{intern}})\left[\text{mS}/\text{c}{\text{m}}^2\right]=g_{\text{AMPA}}$, with N_intern the number of interneurons in the sPING-network. In each case a single pulse stimulation of 300 ms was applied with a total simulation time of 500 ms. The irradiance was set such that the firing rates were equal for both ChR2 models. The study was performed with both a fixed step (10 μs) runge-kutta 4 solver and an ode15s-solver (stiff VSVO-solver, order 1-5, based on numerical differentiation formulas) (Shampine and Reichelt, 1997) with a maximum time step of 100 μs, and a relative and absolute tolerance of 10⁻⁶.

The results shown in this paper are computed with a 3.4 GHz clock rate, quad core system and 8 GB RAM.

2.4. Versatility

The versatility of the proposed model structure is shown with a fit to the MerMAID1 opsin (Oppermann et al., 2019). For more detail on the data set, we refer to the work of Oppermann et al. (2019). The same metrics as aforementioned are used to assess the fit accuracy.

3. Results

To test the feasibility of the proposed double two-state opsin model structure (22OM), it was fit to two data sets. First, we fitted the model to the data set of a ChR2(H134R) opsin reported by Williams et al. (2013), which was collected in a ChR2(H134R)-HEK293 stable cell line (Williams et al., 2013). By the same group already a four state Markov model was fit. This allowed us to analyze the performance of our model in detail. To this end, a comparison of the response to optical stimuli was made in a regular spiking neuron (Pospischil et al., 2008). Moreover, the computational speed was determined for different stimulation paradigms in the former neuron model as well as in the sPING (Börgers and Kopell, 2005) network model with increasing number of transfected neurons. Finally the versatility of the proposed modeling scheme was assessed with a fit to a MerMAID opsin which is an anion-conducting and intensely desensitizing channelrhodopsin.

3.1. The ChR2(H134R) Fit

A 22OM fit of the ChR2(H134R) opsin was obtained by applying the fitting procedure, described in the materials and methods section 2.2, to the experimental data. As Williams et al. (2013) already reported the target features, the first step could be omitted. The absence of differential equations in our fitting procedure allowed for multiple fits to be made, due to the significant reduction of the computational cost. Multiple weight sets, non-linear constraints and combinations of dependency addition of the time constants (product Equation 15 and reciprocal sum Equation 16) were tested. The parameters of the two best fits are shown in Table 1, where RSRS and PP is the fit with a double reciprocal sum and product combination, respectively. Both results were obtained with w_peak = 10, w_ss = 20, w_ratio = 50, w_on = 1000, w_inact = 1000, w_off = 1000, w_recov = 20, and a constraint where O_∞(I, V) > 0.6 for I ≥ 5500 W/m². The weights are chosen as such to level the differences between features to the same order of magnitude. As a result, all features have the same impact in the cost function with a slight preference for the current features. The time constant features are all expressed in seconds, while their values are in the order of milliseconds (except τ_recov), explaining the high weight values. The extra constraint is justified as the current peak already starts to saturate for the highest intensity values, thus clamping the intensity dependence of the open steady-state value above the bending point in the logistics curve.

The models' accuracy according to the four goodness-of-fit metrics Equation 22 are shown in Figure 2. Overall, a positive effect of the final optimization step can be observed. The largest impact is on the time constants, as expected. In the second step of the fitting procedure, the transition rate time-constants (τ_O and τ_R) are approximated with a one on one relationship of the target features (see Equation 9 and 10). These approximations are true in case of high differences in order of magnitude. However, when the differences are smaller some cross correlations exist, for instance τ_R strongly affects τ_on as well, resulting in an underestimation of τ_O. We denote that according to all metrics, the estimation accuracy of τ_on and τ_inact increases, however, at the cost of τ_off. Also, a significant improvement is observed in case of τ_recov. This deviation is due to the fact that the conditions Equation 11 are not fully met. Furthermore, an increased goodness-of-fit of the inactivation time constant can be observed in case of the RSRS vs PP fit. τ_R predominately defines both the inactivation and recovery time constant. In case of the PP fit, a separation of variables is applied where independence is assumed. However, as can be seen in Figures 3F,H, a more clear voltage dependency is present in τ_recov compared to τ_inact. In other words, for low intensities (with high time constants as result) the potential effect is high while the effect is low for high intensities or small time constants. This interdependence is exactly obtained with the reciprocal addition scheme. The same, however less pronounced, can be observed in case of the activation and deactivation time constants (τ_on and τ_off). Consequently, only the RSRS fit is used in further analysis.

FIGURE 2

Figure 2. Normalized goodness-of-fit results of the model fits to the ChR2(H134R) data reported by Williams et al. (2013). Goodness-of-fit according to four metrics: Root-mean-square error (A), Root-mean-square normalized error (B), Root-mean-square weighted error (C) and Root-mean-square of z-score error (D). The compared models are the 4SB model reported by Williams et al. (2013), the 22OM model with twice the product combination of time constants (intermediate fit 22OM: PP-intm and final fit 22OM: PP-final), and the 22OM model with twice the reciprocal addition combination of time constants (intermediate fit 22OM: RSRS-intm and final fit 22OM: RSRS-final). The normalization value is given in the title of each sub-figure. all, τ_all and I_all, are errors where the squared errors of all features, all time constants and all current features are, respectively, added first before taking the root and mean. The depicted legend is valid in all sub-figures.

FIGURE 3

Figure 3. Comparison of model outcomes (4SB and 22OM: RSRS-final) with parameters obtained from experiments. (A), The ChR2(H134R) current during a pulse of 0.5 s (indicated by blue bar) at a voltage clamp of –60 mV; according to the 4SB model (full lines) and 22OM model (dashed dotted lines). The colors indicate the applied intensity and are valid in (A–G). The dotted line and square indicate respectively the experimental current peak and steady-state current at corresponding intensity and potential. (B,D,F), Voltage dependence of respective τ_on, τ_off, and τ_inact across four irradiance levels. (C,E,G), The current-voltage curves of the peak, steady-state and current ratio, respectively. The asterisks with errorbars indicate the experimental mean ± standard deviation. (H), The recovery time constant as function of the membrane potential for three different irradiance levels as depicted in the plot.

Figure 3 shows a detailed comparison of the outcome of our model according to the RSRS fit and the 4SB model, vs. the experimentally determined target features. Overall, it can be observed that the proposed model performs at least as well as the 4SB model. Moreover, all features are well approximated. It can be seen that with the 4SB model, the steady-state value is overestimated in case of negative potentials (Figures 3A,E). However, a better representation is obtained for positive potentials, which explains the lower root-mean-squared normalized error (RMSNE, Figure 2B).

3.2. Neural Response in Regular Spiking Neuron

To analyze the neural response, the strength duration curves (SDC) are determined of the proposed 22OM model with RSRS fit and the 4SB model in a regular spiking neuron, described in Pospischil et al. (2008). First, the Hill-Lapicque model fit is performed on temporal average current (TAC), as described in section 2.3. Very good fits were obtained for both models. The adjusted r² (${\stackrel{-}{\text{R}}}^2$) of TAC vs. PD are 0.9961 and 0.9953 for the 22OM and 4SB model, respectively. The rheobase of the 22OM model (0.49 μA/cm²) is slightly higher than when the 4SB is used (0.47 μA/cm²). Also the chronaxie is higher (47.51 ms vs. 39.45 ms). Consequently, according to the 4SB model for any pulse duration, less charge is injected optogenetically to excite a regular spiking neuron via a ChR2(H134R) opsin. The difference between the models can be attributed to the difference in deactivation time constant (τ_off). This is higher in the 22OM model resulting in a slower closing mechanism and thus increased current injection after the AP. A good cell-type-specific empirical mapping of TAC to irradiance was obtained as well (Equation 25), with ${\stackrel{-}{\text{R}}}^2$ values of 0.9449 (22OM) and 0.9638 (4SB). The parameter values are respectively, a = 8.18, b = 1.26 and c = 1.68, and a = 22.30, b = 1.51 and c = 12.32 in case of the 22OM and 4SB mapping. The lower ${\stackrel{-}{\text{R}}}^2$ of the 22OM mapping resulted also in a slightly lower value of 0.9298 for the irradiance to PD curve while this is 0.9509 in case of the 4SB fit. Based on the mapping parameters and Figure 4, it can be seen that lower intensity level results in higher injected currents when the 22OM model is used. Indeed, extrapolation of the model fit to low intensities results in higher open probabilities than for the 4SB model, hence the difference in irradiance rheobase of 4.90 W/m² vs. 19.01 W/m². Based on the higher peak values for high intensities in case of the 4SB model, one could expect convergence of the irradiance SDCs. However, due to the slow activation kinetics, the peak value is not reached at small pulse durations. Even though the activation time constant is overall higher for the 22OM model (Figure 3B), the bi-exponential current rise due to the extra state variable [τ_ChR2·dp/dt = S0(I) − p, a time-dependent function reflecting the probabilistic, non-instantaneous response of the ChR2-retinal complex to light Williams et al., 2013] in the 4SB model results in a lower current value at the end of the pulse.

FIGURE 4

Figure 4. The strength duration curves (SDC) of the 22OM RSRS and 4SB model in a regular spiking neuron. (A), Irradiance vs. pulse duration with a mapping (dashed line) of the SDC in (B) according to a power series. (B), Temporal average current or average injected current vs. pulse duration. Dashed line represents the Hill-Lapicque model fit. The rheobase and chronaxie are depicted in the figures. The results of the 22OM and 4SB model are in purple and green, respectively.

3.3. Computational Speed

The proposed model in this study contains only two differential equations, which is 50% less in comparison with the 4SB model. Consequently, a reduction of the computational time is expected. Figures 5A–F summarizes the computational speed for different stimulation protocols in a regular spiking neuron. This for fixed irradiances (22OM: 3162 W/m² and 4SB: 1259 W/m²) set to a value that elicit a firing rate of 100 Hz, as described in section 2.3. Figures 5A–D show an overall increase of the computational speed in favor of the 22OM model, with a maximum of 25% for high frequency and duty cycle stimulation. On average the relative difference of the simulation speeds, i.e., simulation speed with 22OM minus simulation speed with 4SB with respect to the latter, is about 20%. Because the simulations were solved using a variable step solver, the difference in firing rate could distort the effective simulation speed, as during an action potential a smaller timestep is selected. Therefore, the relative difference of the simulation speed normalized to the firing rate is depicted as well, with an increase of the gain to 60% as result. The runtime vs. number of transfected neurons is depicted in Figures 5G–I. The simulation outcomes were the same with the variable and fixed step solver, validating the solver settings. Moreover, the firing rate was equal for both opsin models, hence no normalization was necessary. A clear reduction can be observed when the 22OM model is selected instead of the 4SB model, both with a fixed and variable step solver. The time gain by using the proposed model is 15% (5%) in case of 12 neurons and goes up to 40% (15%) and rising when 400 transfected neurons are included with a variable (fixed) step solver.

FIGURE 5

Figure 5. The computational speed of optogenetic neuromodulation in a regular spiking (RS) neuron and sparse Pyramidal-Interneuron-Network-Gamma (sPING). (A–F), Simulation speed, i.e., simulation time/runtime, for different stimulation protocols with varying pulse duration (PD) and pulse repetition frequency (PRF) in a regular spiking neuron, described by Pospischil et al. (2008). (A), The absolute simulation speed with the 22OM-RSRS fit. (B), The simulation speed with the 4SB model. Colorbar is valid for (A) and (B). (C), The relative difference, i.e., (22OM-4SB)/4SB. (D), The effect of the duty cycle on the simulation speed. (E), The difference in firing rate in case of the 22OM model vs. 4SB. (F), The relative difference of simulation speed normalized to the firing rate. (G–I), Runtime of a continuous 300 ms optical pulse in the sparse Pyramidal-Interneuron-Network-Gamma (sPING), with increasing number of transfected neurons. (G), Runtime with a variable step solver. (H), Runtime with a fixed step solver. (I), Relative computation gain, i.e., -(22OM-4SB)/4SB. The used intensities are shown in the titles of (A,B), which give rise to a 100 Hz firing rate (see section 2.3).

3.4. Versatility of the Proposed Model

Finally, we address the versatility of the proposed model and the fitting procedure. Due to the increasing number of possible opsins, it is favorable that their kinetics can be correctly modeled and a fit is easily obtained without preliminary knowledge. To this end, we applied our fitting procedure to experimental data of a MerMAID opsin, which has unlike classical ChR2 a very strong desensitization (Oppermann et al., 2019). Starting from the photocurrent traces, the target features had to be extracted first. Next the parameter space was defined. The rectification function was omitted because this was not observed in the experimental data. Aside from this, the lower bound and initial condition of only the third parameter of R_∞ was altered (Table 1). This straight forward adjustment was made due to the strong desensitization. The weights of the cost function were set to w_peak = 0.04, w_ss = 1, w_ratio = 250, w_on = 10000, w_inact = 10000, w_off = 10000, w_recov = 10, again to level the errors to the same order of magnitude. Because no saturation of the current was observed at high intensity levels a constraint: O_∞(I, V) < 0.5 for I ≤ 4000W/m², was added.

The result of the fit is shown in Figure 6. The parameters of the final and intermediate fit are summarized in Table 1. The model here is with a double product combination of the time constant dependencies. Because the recovery time constant was only determined under one condition, there is no evidence on the interdependence of the variables. This is also supported by the small voltage dependence of the (de)activation time constants. Overall, it can be stated that a good fit is obtained as all kinetics are expressed correctly. Only, the deactivation time constant seems to be underestimated. This is a consequence of the constraint in Equation 21, which ensures a current decay back to baseline after optical stimulation. Due to its strong desensitization, R_∞ has to be small, thus inducing an upper limit on τ_O(0, V), which defines the deactivation time constant. The trade off is justified due to the higher uncertainty of the deactivation time constant (see Figure 6D). Moreover, the overall effect is expected to be low as can be seen in the right inset of Figure 6A.

FIGURE 6

Figure 6. Comparison of the 22OM-Mermaid (final fit) model outcomes and experimental data. (A), In gray, the photocurrent of a voltage clamp experiment during a 0.5 s continuous illumination with an intensity of 3,734 W/m² (indicated with blue bar on top) (Oppermann et al., 2019); In red, the corresponding model outcome. Left inset is a zoom of the current peak (0.045–0.075 s, indicated with black bar). Right inset is a zoom of the current deactivation (0.45–0.7 s, indicated with a blue square). (B–G), The voltage dependence of the target features (τ_on, I_peak, τ_off, I_ss, τ_inact, and I_ratio) at an irradiance of 3,734.4 W/m² is shown in blue. The light dependence at a holding potential of –60 mV is depicted in red. (H), Ratio of the peak currents in response to a two-pulse stimulation protocol at –60 mV and 3,734 W/m² as function of the inter-pulse interval. The recovery time (the interval time necessary to have a ratio of 63%), is indicated with a black arrow.

4. Discussion

The proposed double two-state model structure for the modeling of opsins appears to be a good alternative to the computationally more expensive four state Markov, non-instantaneous models. All features are represented, with even some improved fit accuracy in comparison with a four state Markov variant. Furthermore, with the proposed fitting procedure, we were able to fit two opsins, ChR2(H134R) and MerMAID. Although the prominent difference of the mutants kinetics, the fitting procedure allowed us to get these fits with only minor adjustments of the parameter space and constraints. Therefore, creating the possibility for autonomous model fitting based on photocurrent traces. Moreover, a good fit is obtained within an acceptable time frame, due to the absence of differential equations in the fitting procedure, which is not achievable in case of a four state Markov model. The intermediate fit is obtained within 3 h, while the final fit always flagged the time limit of 24 h. Increasing the limit improves the fit accuracy but only small changes were observed. Fine tuning of the optimization settings, such as number of particles or tolerances, could reduce the training error even more. However, this is out of the scope of this study.

The proposed model is an empirical model. The fit is performed on a limited dataset thus extrapolation should be treated with care. This is clear from the neural response results in section 3.2. Although both the 4SB and our model were fit to the same experimental data, a clear discrepancy between the fitted rheobase is observed [4.90 W/m² (22OM) vs. 19.01 W/m² (4SB)]. Unlike the chronaxie where the difference can be attributed to the model's structure, the difference in rheobase is due to the discrepancy between opening rates after extrapolation to low intensities, attributed to the fit and intensity dependence chosen in each model. More experiments are required in order to validate this.

The dependencies chosen here are all, except the rectification, sigmoidal. Therefore, they are all bounded and monotonic. This is in accordance with a channel's behavior, i.e., increased and faster opening at higher intensities but limited to an open probability of one. We opted for a biphasic logistics function for τ_R(I) modeling. This is in agreement with the hypothesis of the necessity of two light dependent rates [(R → S) and (S → R), see section 2.1] and the second and third photochemical pathways described by Kuhne et al. (2019) (Figure 1). Other functions were tested, e.g., weibull or asymmetric logistics with double intensity dependence, however no improvement was observed. Initially, separation of variables was assumed to suffice due to the lack of experimental evidence of complex channel interdependence of both irradiance and potential of each feature separately. However, due to the models structure, τ_on and τ_off share the same voltage dependency, as well as τ_inact an τ_recov. The voltage dependence of τ_recov and τ_off was clearly more pronounced in the experimental data of the ChR2(H134R) mutant. Therefore, the reciprocal addition Equation 16 was tested as alternative, resulting in an improved fit accuracy. However, this only scales down the voltage dependent effect on τ_on and τ_inact while the same relationship is maintained. The necessity of more complex relationships could be investigated in future work as well as the need for voltage dependence of the rate functions steady-state values (O_∞ and R_∞), which was omitted in this study.

Currently the model incorporates voltage and irradiance dependence. Studies have however shown the importance of pH on the channel kinetics in many opsins. Furthermore, ion concentrations have an impact on the reversal potentials and current rectification (Berndt et al., 2010; Stehfest and Hegemann, 2010). Schneider et al. (2013), postulated a model based on the kinetics of multiple ion species interacting with the channel, with an improved representation of the current rectification (Schneider et al., 2013; Williams et al., 2013). While the photocurrent properties are unaffected by pH-changes, the MerMAID photocurrent is strongly dependent on the Cl⁻ concentrations. The fit performed here was on experimental data recorded with an extracellular Cl⁻ concentration of 150 mM and intracellular Cl⁻ of 120 mM, explaining the depolarizing currents (negative sign in Figure 6) as an anion conducting channel. By changing the extracellular concentration to 10 mM, the channel's reversal potential is shifted to the reversal potential of Cl⁻. (The concentrations are exchanged with respect to a conventional neuron, where the typical intracellular and extracellular concentrations are 10 and 120 mM, respectively. This explains the experimentally measured depolarizing currents (negative sign), while one would expect hyperpolarizing currents (positive sign) from a Cl⁻ conducting channel.) Evidence of the Cl⁻ effect on channel kinetics is still absent but further experiments are needed (Oppermann et al., 2019). Consequently, the model fit shown here can be used in computational studies but the reversal potential should be adjusted accordingly.

With the current model structure, the model responds instantaneously to light (see left inset Figure 6A). With the 4SB model this is circumvented by adding a extra state variable with a time constant of 1.5 ms. It is clear that for long (PD >> τ_on) continuous pulses its effect is negligible, as activation is dominated by the activation time constant. However, with short bursts or pulses, this non-instantaneous activation becomes prominent as observed in section 3.2. In future work, it could therefore be interesting to incorporate this non-instantaneous response. This could probably be obtained by adding an extra state variable, as performed with the 4SB model, however at the cost of the computational speed. Another possibility is to raise the open state, O(t), to a higher power, smoothing the transition but without irradiance control. Modification of the model's structure could be circumvented by gradually increasing the intensity, instead of applying a rectangular pulse.

5. Conclusion

To facilitate computational studies in the field of optogenetics, we proposed a double two-state opsin model structure as alternative to the conventional three and four state Markov models. In the proposed model, the second state pair represents the conductance regulation due the dark-light adaptation. With this model type, a reduction in complexity is obtained resulting in only two differential equations compared to four in case of the preferred, non-instantaneous four state Markov models used for opsin modeling. With the provided fitting procedure, autonomous model fits of two distinctive opsins (ChR2(H134R) and MerMAID) were obtained. Both model fits were performed within an acceptable time frame thanks to the absence of differential equations and parameter space reduction associated with the multi step approach. Moreover, both models are able to represent the experimental data with great accuracy. Due to the model's structure, there is, however, an instantaneous response to light, overestimating the injected current at very short pulses (< τ_on). Furthermore, pH and ion concentration dependence are not incorporated. In its current state with only two differential equations, the computational speed is increased up to 25% in a regular spiking neuron and up to 40% in a network of 400 transfected neurons.

Data Availability Statement

The original contributions presented in the study are publicly available. Requests to access these datasets should be directed to Ruben.Schoeters@UGent.be.

Author Contributions

RS, TT, LM, WJ, RR, and ET contributed to the study design. RS performed all simulations. RS, TT, and ET prepared the manuscript. All authors read and approved the final manuscript.

Funding

RS is a Ph.D. Fellow of the FWO-V (FR) (Research Foundation Flanders, Belgium; grant no. 11E4520N). TT is a PhD Fellow of the FWO-V (SB) (Research Foundation Flanders, Belgium; grant no. 1S00417N).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Johannes Oppermann for sharing his experimental data on MerMAIDs. Exploration of the parameter space and dependencies in this work was carried out using the Supercomputer Infrastructure (STEVIN) at Ghent University, funded by Ghent University, the Flemish Supercomputer Center (VSC), the Hercules Foundation and the Flemish Government department EWI.

References

Abilez, O. J., Wong, J., Prakash, R., Deisseroth, K., Zarins, C. K., and Kuhl, E. (2011). Multiscale computational models for optogenetic control of cardiac function. Biophys. J. 101, 1326–1334. doi: 10.1016/j.bpj.2011.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Aravanis, A. M., Wang, L. P., Zhang, F., Meltzer, L. A., Mogri, M. Z., Schneider, M. B., et al. (2007). An optical neural interface: in vivo control of rodent motor cortex with integrated fiberoptic and optogenetic technology. J. Neural Eng. 4:S143–S156. doi: 10.1088/1741-2560/4/3/S02

PubMed Abstract | CrossRef Full Text | Google Scholar

Azimihashemi, N., Erbguth, K., Vogt, A., Riemensperger, T., Rauch, E., Woodmansee, D., et al. (2014). Synthetic retinal analogues modify the spectral and kinetic characteristics of microbial rhodopsin optogenetic tools. Nat. Commun. 5, 1–12. doi: 10.1038/ncomms6810

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, A., Prigge, M., Gradmann, D., and Hegemann, P. (2010). Two open states with progressive proton selectivities in the branched channelrhodopsin-2 photocycle. Biophys. J. 98, 753–761. doi: 10.1016/j.bpj.2009.10.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Börgers, C., and Kopell, N. (2005). Effects of noisy drive on rhythms in networks of excitatory and inhibitory neurons. Neural Comput. 17, 57–608. doi: 10.1162/0899766053019908

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruun, S., Stoeppler, D., Keidel, A., Kuhlmann, U., Luck, M., Diehl, A., et al. (2015). Light-dark adaptation of channelrhodopsin involves Photoconversion between the all-trans and 13-cis retinal isomers. Biochemistry 54:5400. doi: 10.1021/acs.biochem.5b00597

CrossRef Full Text | Google Scholar

Carrette, S., Boon, P., Sprengers, M., Raedt, R., and Vonck, K. (2015). Responsive neurostimulation in epilepsy. Exp. Rev. Neurotherapeut. 15, 1445–1454. doi: 10.1586/14737175.2015.1113875

CrossRef Full Text | Google Scholar

Chen, S. (2018). Constrained Particle Swarm Optimization. MATLAB File Exchange.

Chen, Y., Xiong, M., and Zhang, S. C. (2015). Illuminating Parkinson's therapy with optogenetics. Nat. Biotechnol. 33, 149–150. doi: 10.1038/nbt.3140

CrossRef Full Text | Google Scholar

Deisseroth, K., and Hegemann, P. (2017). The form and function of channelrhodopsin. Science 357:eaan5544. doi: 10.1126/science.aan5544

CrossRef Full Text | Google Scholar

Ernst, O. P., Sánchez Murcia, P. A., Daldrop, P., Tsunoda, S. P., Kateriya, S., and Hegemann, P. (2008). Photoactivation of channelrhodopsin. J. Biol. Chem. 283, 1637–1643. doi: 10.1074/jbc.M708039200

CrossRef Full Text | Google Scholar

Gerits, A., and Vanduffel, W. (2013). Optogenetics in primates: a shining future? Trends Genet. 29, 403–411. doi: 10.1016/j.tig.2013.03.004

CrossRef Full Text | Google Scholar

Grossman, N., Simiaki, V., Martinet, C., Toumazou, C., Schultz, S. R., and Nikolic, K. (2013). The spatial pattern of light determines the kinetics and modulates backpropagation of optogenetic action potentials. J. Comput. Neurosci. 34, 477–488. doi: 10.1007/s10827-012-0431-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunaydin, L. A., Yizhar, O., Berndt, A., Sohal, V. S., Deisseroth, K., and Hegemann, P. (2010). Ultrafast optogenetic control. Nat. Neurosci. 13, 387–392. doi: 10.1038/nn.2495

CrossRef Full Text | Google Scholar

Hassan, R., Cohanim, B., de Weck, O., and Venter, G. (2005). “A Comparison of Particle Swarm Optimization and the Genetic Algorithm,” 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference (Austin, TX), 1–13.

PubMed Abstract

Helwig, S., and Wanka, R. (2007). “Particle swarm optimization in high-dimensional bounded search spaces,” in Proceedings of the 2007 IEEE Swarm Intelligence Symposium, SIS 2007 (Honolulu, HI), 198–205.

Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). Propagation of electrical signals along giant nerve fibres. Proc. R. Soc. Lond. B Biol. Sci. 140, 177–183. doi: 10.1098/rspb.1952.0054

CrossRef Full Text | Google Scholar

Klapoetke, N. C., Murata, Y., Kim, S. S., Pulver, S. R., Birdsey-Benson, A., Cho, Y. K., et al. (2014). Independent optical excitation of distinct neural populations. Nat. Methods 11, 338–346. doi: 10.1038/nmeth.2836

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhne, J., Vierock, J., Tennigkeit, S. A., Dreier, M. A., Wietek, J., Petersen, D., et al. (2019). Unifying photocycle model for light adaptation and temporal evolution of cation conductance in channelrhodopsin-2. Proc. Natl. Acad. Sci. U.S.A. 116, 9380–9389. doi: 10.1073/pnas.1818707116

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagel, G., Ollig, D., Fuhrmann, M., Kateriya, S., Musti, A. M., Bamberg, E., et al. (2002). Channelrhodopsin-1: a light-gated proton channel in green algae. Science 296, 2395–2398. doi: 10.1126/science.1072068

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagel, G., Szellas, T., Huhn, W., Kateriya, S., Adeishvili, N., Berthold, P., et al. (2003). Channelrhodopsin-2, a directly light-gated cation-selective membrane channel. Proc. Natl. Acad. Sci. U.S.A. 100, 13940–13945. doi: 10.1073/pnas.1936192100

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikolic, K., Grossman, N., Grubb, M. S., Burrone, J., Toumazou, C., and Degenaar, P. (2009). Photocycles of channelrhodopsin-2. Photochem. Photobiol. 85, 400–411. doi: 10.1111/j.1751-1097.2008.00460.x

CrossRef Full Text | Google Scholar

Noble, D., and Stein, R. B. (1966). The threshold conditions for initiation of action potentials by excitable cells. J. Physiol. 187, 129–162. doi: 10.1113/jphysiol.1966.sp008079

CrossRef Full Text | Google Scholar

Oppermann, J., Fischer, P., Silapetere, A., Liepe, B., Rodriguez-Rozada, S., Flores-Uribe, J., et al. (2019). MerMAIDs: a family of metagenomically discovered marine anion-conducting and intensely desensitizing channelrhodopsins. Nat. Commun. 10, 1–13. doi: 10.1038/s41467-019-11322-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., et al. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol. Cybernet. 99, 427–441. doi: 10.1007/s00422-008-0263-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritter, E., Piwowarski, P., Hegemann, P., and Bartl, F. J. (2013). Light-dark adaptation of channelrhodopsin C128T mutant. J. Biol. Chem. 288, 10451–10458. doi: 10.1074/jbc.M112.446427

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, F., Gradmann, D., and Hegemann, P. (2013). Ion selectivity and competition in channelrhodopsins. Biophys. J. 105, 91–100. doi: 10.1016/j.bpj.2013.05.042

CrossRef Full Text | Google Scholar

Schoeters, R., Tarnaud, T., Joseph, W., Martens, L., Raedt, R., and Tanghe, E. (2019). “Channelrhodopsin-2 model with improved computational efficiency,” in 28th Annual Computational Neuroscience Meeting: CNS*2019, Vol. 20 (Barcelona: BMC Neuroscience), 194.

Google Scholar

Shampine, L. F., and Reichelt, M. W. (1997). The matlab ode suite. SIAM J. Sci. Comput. 18, 1–22. doi: 10.1137/S1064827594276424

CrossRef Full Text | Google Scholar

Sherfey, J. S., Soplata, A. E., Ardid, S., Roberts, E. A., Stanley, D. A., Pittman-Polletta, B. R., et al. (2018). DynaSim: a MATLAB toolbox for neural modeling and simulation. Front. Neuroinform. 12:10. doi: 10.3389/fninf.2018.00010

PubMed Abstract | CrossRef Full Text | Google Scholar

Stehfest, K., and Hegemann, P. (2010). Evolution of the channelrhodopsin photocycle model. ChemPhysChem 11, 1120–1126. doi: 10.1002/cphc.200900980

CrossRef Full Text | Google Scholar

Tønnesen, J., and Kokaia, M. (2017). Epilepsy and optogenetics: can seizures be controlled by light? Clin. Sci. 131, 1605–1616. doi: 10.1042/CS20160492

PubMed Abstract | CrossRef Full Text | Google Scholar

Volkov, O., Kovalev, K., Polovinkin, V., Borshchevskiy, V., Bamann, C., Astashkin, R., et al. (2017). Structural insights into ion conduction by channelrhodopsin 2. Science 358:eaan8862. doi: 10.1126/science.aan8862

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, J. C., and Entcheva, E. (2015). Optogenetic versus electrical stimulation of human cardiomyocytes: Modeling insights. Biophys. J. 108, 1934–1945. doi: 10.1016/j.bpj.2015.03.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, J. C., Xu, J., Lu, Z., Klimas, A., Chen, X., Ambrosi, C. M., et al. (2013). Computational optogenetics: empirically-derived voltage- and light-sensitive channelrhodopsin-2 model. PLoS Comput. Biol. 9:e1003220. doi: 10.1371/journal.pcbi.1003220

PubMed Abstract | CrossRef Full Text | Google Scholar