Double Two-State Opsin Model With Autonomous Parameter Inference

Optogenetics has a lot of potential to become an effective neuromodulative therapy for clinical applications. Selecting the correct opsin is crucial to have an optimal optogenetic tool. With computational modeling, the neuronal response to the current dynamics of an opsin can be extensively and systematically tested. Unlike electrical stimulation where the effect is directly defined by the applied field, the stimulation in optogenetics is indirect, depending on the selected opsin's non-linear kinetics. With the continuous expansion of opsin possibilities, computational studies are difficult due to the need for an accurate model of the selected opsin first. To this end, we propose a double two-state opsin model as alternative to the conventional three and four state Markov models used for opsin modeling. Furthermore, we provide a fitting procedure, which allows for autonomous model fitting starting from a vast parameter space. With this procedure, we successfully fitted two distinctive opsins (ChR2(H134R) and MerMAID). Both models are able to represent the experimental data with great accuracy and were obtained within an acceptable time frame. This is due to the absence of differential equations in the fitting procedure, with an enormous reduction in computational cost as result. The performance of the proposed model with a fit to ChR2(H134R) was tested, by comparing the neural response in a regular spiking neuron to the response obtained with the non-instantaneous, four state Markov model (4SB), derived by Williams et al. (2013). Finally, a computational speed gain was observed with the proposed model in a regular spiking and sparse Pyramidal-Interneuron-Network-Gamma (sPING) network simulation with respect to the 4SB-model, due to the former having two differential equations less. Consequently, the proposed model allows for computationally efficient optogenetic neurostimulation and with the proposed fitting procedure will be valuable for further research in the field of optogenetics.


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). Postillumination, there is a bi-exponential decay back to baseline, rendering the channel in an apparent dark-adapted state (DA app ).
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. 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).
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

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.

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: 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: with the Heaviside function, O 0 and R 0 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.
Furthermore, following conditions need to be met for the relationship to hold true: 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 2 . The units of the conductance depend on the experimental data of each opsin, i.e., mS/cm 2 and µS in case of the ChR2(H134R) and MerMAID fit, respectively.

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 multistart algorithm with ten starting points was used to ensure finding of the global solution. The variable and function tolerance were set to 10 −12 . 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.
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: Here, y x and t x,I i ,V i are, respectively the model output and target value at stimulation values (I,V), with y peak = max t (|i on ChR2 (t, I, V)|), y ss = i on ChR2 (t off ) and y ratio = y ss /y peak . i on ChR2 (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 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: τ on (I, V) 2 , τ off (I, V) 2 , τ inact (I, V) 2 and τ 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 −4 s, while for the two-pulse experiments a step of 1 ms is used.

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): where w x equals 1, 1/t x,I i ,V i , or 1/σ x,I i ,V i in case of RMSE, RMSNE or RMSZE, respectively. y x (I i , V i ), t x,I i ,V i and σ x,I i ,V i 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,  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).
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: 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 ChR2,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) ode113solver (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 −3 and 10 −6 , 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 GABAa = 2/(N intern )[mS/cm 2 ] = g 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 −6 .
The results shown in this paper are computed with a 3.4 GHz clock rate, quad core system and 8 GB RAM.

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.

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.

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,  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 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).

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 2 (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 2 ) is slightly higher than when the 4SB is used (0.47 µA/cm 2 ). 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), withR 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 lowerR 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 2 vs. 19.01 W/m 2 . 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.

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 2 and 4SB: 1259 W/m 2 ) 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.

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 2 , 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.

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, noninstantaneous 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  (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 2 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 2 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.
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 2 (22OM) vs. 19.01 W/m 2 (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 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.

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.