Home | About | Journals | Submit | Contact Us | Français |

**|**Front Neuroinform**|**v.10; 2016**|**PMC4834562

Formats

Article sections

Authors

Related links

Front Neuroinform. 2016; 10: 8.

Published online 2016 March 11. doi: 10.3389/fninf.2016.00008

PMCID: PMC4834562

Edited by: Andrew P. Davison, Centre National de la Recherche Scientifique, France

Reviewed by: Padraig Gleeson, University College London, UK; Emilia Entcheva, Stony Brook University, USA

Received 2015 November 4; Accepted 2016 February 19.

Copyright © 2016 Evans, Jarvis, Schultz and Nikolic.

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Optogenetics has become a key tool for understanding the function of neural circuits and controlling their behavior. An array of directly light driven opsins have been genetically isolated from several families of organisms, with a wide range of temporal and spectral properties. In order to characterize, understand and apply these opsins, we present an integrated suite of open-source, multi-scale computational tools called PyRhO. The purpose of developing PyRhO is three-fold: (i) to characterize new (and existing) opsins by automatically fitting a minimal set of experimental data to three-, four-, or six-state kinetic models, (ii) to simulate these models at the channel, neuron and network levels, and (iii) provide functional insights through model selection and virtual experiments *in silico*. The module is written in Python with an additional IPython/Jupyter notebook based GUI, allowing models to be fit, simulations to be run and results to be shared through simply interacting with a webpage. The seamless integration of model fitting algorithms with simulation environments (including NEURON and Brian2) for these virtual opsins will enable neuroscientists to gain a comprehensive understanding of their behavior and rapidly identify the most suitable variant for application in a particular biological system. This process may thereby guide not only experimental design and opsin choice but also alterations of the opsin genetic code in a neuro-engineering feed-back loop. In this way, we expect PyRhO will help to significantly advance optogenetics as a tool for transforming biological sciences.

Optogenetics is a biotechnology which renders excitable cells light-sensitive by inserting genes which, upon expression, create light-activated ion channels known originally as rhodopsins (Nagel et al., 2003; Boyden et al., 2005). Over the last 10 years optogenetics has found widespread application, initially in neuroscience (Zhang et al., 2006; Adamantidis et al., 2007; Arenkiel et al., 2007; Han and Boyden, 2007; Yizhar et al., 2011), but increasingly also in more distal areas of physiology such as cardiac science (Arrenberg et al., 2010; Boyle et al., 2013), intracellular signaling (Airan et al., 2009), and gene transcription (Konermann et al., 2013). Applications to date have included control of motor cortex (Aravanis et al., 2007), cortical circuit mapping (Wang et al., 2007; Zhang et al., 2007; Ayling et al., 2009; Petreanu et al., 2009; Klapoetke et al., 2014), optoelectronic neuroprosthetic devices (for example retinal Lagali et al., 2008; Degenaar et al., 2009; Busskamp and Roska, 2011 or cochlear prostheses, Hernandez et al., 2014), regulation of the symptoms of neurodegenerative disorders (e.g., Parkinson's Gradinaru et al., 2009), closed loop control of epileptic seizures, peripheral nerve stimulation (Arlow et al., 2013), and novel cardiac pacemaker technology (Bruegmann and Sasse, 2015), to name just a few.

In an effort to develop more effective and tailored opsins, hybrids and genetic mutants are continually being created (Berndt et al., 2008; Lin et al., 2009; Hegemann and Moglich, 2011; AzimiHashemi et al., 2014). Experimentally characterizing these new variants is a lengthy process requiring substantial effort before they can be harnessed to address questions in neuroscience (Chow et al., 2010; Gunaydin et al., 2010; Lin, 2011; Chuong et al., 2014). The problem is further compounded when considering the number of combinations between opsins (with total variations now in the hundreds, Zhang et al., 2011) and target cell types. Experimentally testing each combination of opsin and target cell type of interest is practically impossible, effectively limiting the use of optogenetics as a tool.

Theoretical understanding of the underlying mechanisms of optogenetics has developed over the past 10 years (Hegemann et al., 2005; Feldbauer et al., 2009), which has led to a deeper understanding of the biophysical mechanisms of the photosensitization agents which form the foundations of optogenetics (Hegemann et al., 2005; Bamann et al., 2008; Ernst et al., 2008; Nikolic et al., 2009; Stehfest and Hegemann, 2010). Furthermore, the design and engineering of optogenetic devices must start with models of the underlying molecular mechanisms of opsin behavior in cells (Gradinaru et al., 2007; Nikolic et al., 2007; Shoham and Deisseroth, 2010; Foutz et al., 2012; Williams et al., 2013). Computational modeling is thus core to understanding how light induced ionic transport across cell membranes can be tailored for different applications: from probing cellular physiology to creating new treatments for neurological and psychiatric illnesses.

The quest to both expand and refine optogenetics as an effective tool for neuroscience and other areas of physiology requires multiple levels of analysis: from molecular modeling through kinetic models and even network level models. To aid in this effort we propose *PyRhO*; an integrated suite of open-source, multi-scale computational tools to characterize opsins, then rapidly develop and conduct virtual experiments with them *in silico*.

PyRhO offers several integrated computational tools for analysing and experimenting with (rhod)opsins in a virtual environment:

- The first tool will automatically fit a choice of models to experimental data, extracting the parameters that describe the functional dynamics of the opsins.
- The second tool can then take these extracted parameters (or alternatively use default values) and simulate a wide range of experimental protocols to reproduce the photo-response of the opsin of interest. These protocols are typically voltage-clamp experiments and include common engineering inputs such as steps, ramps, and chirps, along with more tailored protocols such as pairs of increasingly spaced pulses for evaluating the recovery process.
- These models and protocols can be run on several simulation platforms spanning multiple scales (to model isolated opsins or transfected neurons) including:
- Pure Python for simple channel-level voltage clamp experiments;
- NEURON for morphologically detailed models of optogenetically transfected neurons;
- Brian2 for simulating whole networks with transfected groups of neurons.

- A Graphical User Interface (GUI) for easy navigation through all tools, running of virtual experiments and sharing of results.

In this way, PyRhO allows the investigator to simulate opsin dynamics on multiple scales from sub-cellular channels, to individual neurons and finally the dynamics of whole networks. This will help to elucidate the link between the biophysics of opsins and the functional implications of their use in a particular biological system.

The tools are written in Python due to its rapidly growing popularity across the sciences, readability, modularity and large array of open-source modules (Muller et al., 2015). An accompanying GUI running in IPython/Jupyter (Pérez and Granger, 2007) has also been developed to facilitate more interactive exploration of the models for both experimental and pedagogic purposes, requiring virtually no programming experience. In addition to controlling the fitting routines, the GUI also exposes the integrated simulators (e.g., NEURON). Furthermore, this self-logging, notebook-based approach has been identified as a particularly promising medium for sharing models and reproducing results in computational neuroscience (Topalidou et al., 2015).

Simulations based on these virtual opsins will enable neuroscientists to gain insight into their behavior and rapidly identify the most suitable variant for application in a particular biological system, not only guiding choice, but also opsin development. Understanding gained from biologically realistic simulations may provide ideas of how to alter the opsin's genetic code to generate new mutants. These new variations can then be characterized and simulated within PyRhO to determine their suitability for a particular application.

Here, we describe the structure of PyRhO and demonstrate a sample of its capabilities, illustrated through snippets of code and its GUI. We demonstrate the use of PyRhO in fitting models to Channelrhodopsin-2 (ChR2) data and present results for typical illumination strategies and experimental protocols designed to tease apart the effects of key model parameters. We finish with a discussion of the main benefits of using PyRhO, its limitations to date and planned future developments to extend its capabilities.

PyRhO is written as a Python module and released as an open-source project under the revised BSD license. Download and installation instructions can be found with the code at PyRhO's GitHub repository (https://github.com/ProjectPyRhO/PyRhO), along with example notebooks and a link to the project's website containing further information. A virtual machine with all dependencies installed and examples ready to run is also available, such that the GUI can be used with a minimum of set-up and virtually no programming experience.

The module is comprised of several integrated components for fitting model parameters to experimental data and for simulating the models at multiple scales. Fitting data is an optional step since PyRhO is initialized with default parameters, allowing the user to immediately experiment with simulating the three types of opsin models in order to better understand their dynamics. If the required data are provided to the fitting algorithms however, the parameterized models may be run through the stimulation protocols to efficiently characterize the opsins *in silico*, or determine their suitability for a particular application based upon their dynamics.

PyRhO is implemented as a Python package called `pyrho` which builds upon popular scientific Python modules including `scipy`, `numpy`, `matplotlib`, and `lmfit`. Additionally, if optogenetic simulations in detailed morphological models of individual (or a few) neurons are required, `NMODL` files (Hines and Carnevale, 2000) are provided for use with NEURON (Hines et al., 2009). Similarly, for network-level simulations PyRhO has been integrated with the Brian simulator (Goodman and Brette, 2008, 2009) and includes model descriptions suitable for use with Brian2.

The simulation architecture is designed around three layers of abstraction: models, protocols and simulators. These layers are illustrated in the work-flow schematic of Figure Figure11 along with the other major components of PyRhO. Each layer contains families of classes to create a uniform interface for each subclass, for example, the differences in setting the light-dependent transition rates of the three models are shielded from the user by endowing each opsin model subclass with the method `setLight()`. A similar approach is taken with the other layers providing a common set of member variables and methods, making usage consistent and providing a framework for future development of new subclasses (i.e., additional kinetic models, stimulation protocols, and simulation platforms).

A detailed understanding of how the channel is gated and the ions are conducted is still lacking, although some recent studies have elucidated important aspects of the pore formation and ionic transport (Feldbauer et al., 2009; Kuhne et al., 2014) We assume that all light-sensitive ion channel currents (*I*) can be expressed in the classic form:

$$\begin{array}{l}I=g\xb7(v-E)\text{},\end{array}$$

(1)

where *g* is the channel conductance, *v* the membrane voltage and *E* is the reversal potential for the specific opsin type. Generally speaking the ionic conductance is a complex function of light flux (ϕ(*t*)), wavelength (λ), and the opsin's photocycle, membrane voltage, temperature (*T*), and intracellular and extracellular pH (Gradmann et al., 2011). We use a simplified empirical form for the channel conductance, introduced by Hodgkin and Huxley, expressing it as a product of a constant (*g*_{0}, in our case this is maximum conductance at *v* = −70*mV*), and a numerical coefficient (*f* > 0):

$$\begin{array}{l}g=g(\varphi ,\lambda ,v,T,pH,t)={g}_{0}\xb7f(\varphi ,\lambda ,v,T,pH,t)\text{},\end{array}$$

(2)

In this version of PyRhO we have implemented the photocycle and membrane voltage dependencies and assumed that these two contributions can be separated:

$$\begin{array}{l}g={g}_{0}\xb7{f}_{\text{\varphi}}(\varphi ,t)\xb7{f}_{\text{v}}(v)\text{}.\end{array}$$

(3)

These two dependences are considered to be the most relevant for physiological electrolyte conditions, when temperature and pH are considered to be fixed. Other dependencies will be implemented in the next version of PyRhO.

At the core of PyRhO are three functional Markov models of opsin kinetics, namely the three-, four- (Nikolic et al., 2009), and six-state (Grossman et al., 2013) models. We note that very similar models have been investigated in several other studies (Gradmann et al., 2002; Nagel et al., 2003; Hegemann et al., 2005; Ishizuka et al., 2006; Bamann et al., 2008; Ernst et al., 2008; Foutz et al., 2012; Williams et al., 2013) but used our earlier models as a starting point as we have since extended them and unified their notation. These models vary in complexity providing a range in the trade-off between biological accuracy and computational efficiency to choose from. The key features of these models, including an outline of their strengths and weaknesses, are summarized in Table Table11 with accompanying illustrations in Figure Figure2.2. Since their original formulation, the models have been extended to encompass additional parameter dependencies, better fit the experimental data and use a consistent notation, with the full model descriptions given in Table Table2.2. An analytic solution for the three-state model was also calculated and is included in the Appendix.

Both four- and six-state models assume that there are two open states (*O*_{1} and *O*_{2}, see Figure Figure2),2), with channel conductances *g*_{O1} and *g*_{O2}, respectively. The Photocycle factor in Equation (1) has the form:

$$\begin{array}{l}{f}_{\text{\varphi}}(\varphi )={O}_{1}+\gamma {O}_{2}\text{},\end{array}$$

(4)

where *O*_{1} and *O*_{2} are the fractions of opsins in two open states in the interval [0, 1], and γ = *g*_{O2} / *g*_{O1}. In contrast, the three-state model assumes only one open state (*O*) making the photocycle factor simply: *f*_{ϕ}(ϕ) = *O*.

Here, we assume that the membrane voltage affects only the ion-channel conductance but not the channel kinetics. By investigating experimental results for Channelrhodopsin-2 steady-state current vs. clamped voltage, the *I*–*V* curve shows inwardly rectifying behavior (Bamberg et al., 2008). We have previously found that an exponential function gives a good fit for this dependency (Grossman et al., 2011), therefore for the voltage factor in Equation (1) we adopt the form:

$$\begin{array}{l}{f}_{\text{v}}(v)=\frac{{v}_{1}}{v-E}\xb7\left(1-{e}^{-\frac{v-E}{{v}_{0}}}\right)\text{},\end{array}$$

(5)

where *v*_{0} and *v*_{1} are fitting parameters, along with *E*, the channel's reversal potential. The exponential dependence on *v* transforms into a linear dependence for large values of *v*_{0} which cause the exponent to be small and the expression in Equation (5) reduces to *f*_{v}(*v*) ≈ *v*_{1} / *v*_{0} = *const*, i.e., no direct dependence on membrane voltage, which may be a more appropriate form for some opsins. The expression given by Equation (5) therefore generalizes to both cases for appropriate choices of the parameters *v*_{0} and *v*_{1}.

Furthermore, since the voltage dependence factor is defined to be equal to 1 at −70 mV (*f*_{v}(−70 mV): = 1), the value of v_{1} is related to the other parameters through the following equation:

$$\begin{array}{l}{v}_{1}=\frac{70+E}{{e}^{\frac{70+E}{{v}_{0}}}-1}\end{array}$$

(6)

This relationship is used as a constraint in the fitting procedures described below.

PyRhO incorporates novel fitting algorithms with which each of the opsin models may be parameterized when given an appropriate set of data. The fitting algorithms use the `lmfit` module (Newville et al., 2014) which in addition to providing access to a variety of optimization algorithms, enables numerical bounds to be placed on parameter values as well as algebraic constraints. Parameters may also be manually fixed if, for example, they have been directly measured experimentally. Once the algorithm has finished, a `Parameters` object is returned, populated with the extracted parameter values which may then be saved or used as the basis for simulations. Plots of photocurrents simulated with these derived parameters are drawn over the experimental data (with residual error) to allow for visual comparison of the resultant fits.

In order to characterize each model, a set of voltage-clamped photocurrents are required, ideally collected from HEK cells to eliminate the confounding effects of other ion channels which may be present in neurons. To capture all currently modeled variable dependencies, data from three stimulation protocols are necessary, listed below by the model properties which they reveal. In the event of scarce data or uncharacterized variables which the user does not intend to vary, we describe the minimum set of data for the fitting procedure below and discuss the implications for the resultant model.

**Voltage dependence**: {*E, v*_{0},*v*_{1}}. Long light pulse (fixed flux) to steady-state, vary voltage clamp potential (e.g., in steps of 30 mV:*V*_{clamp}= {−100, −70, −40, −10, 20, 50, 80} mV,*n*≥ 5). Voltage clamp values should not be too close to*E*as this may cause distortions in the fitting algorithms. The software will automatically find the plateau values*I*_{ss}, plot*I*_{ss}vs.*V*_{clamp}, and find the fitting parameters for the function*f*_{v}(*v*) given by Equation (5). An example is shown in Figure Figure33.**Recovery rate**: {*G*_{r0}}. Two long light pulses with varying inter-pulse-interval (IPI), Voltage clamp: −70*mV*. Light on (first pulse)—light off for e.g.,*t*_{IPI}= {0.5, 1, 2.5, 5, 10}s—light on (second pulse). The software will automatically find the peak values for each recording, align the data to the end of the first pulse and fit an appropriate exponential curve of the form ${I}_{\text{peak}}(t)={I}_{peak0}-a\xb7{e}^{-{G}_{r0}\xb7{t}_{\text{IPI}}}$. We note here that this expression is strictly speaking correct only when both*O*_{1}and*O*_{2}states are empty. Consequently*a*is left as a free fitting parameter and very short values for*t*_{IPI}should be avoided to prevent the distortions caused by the faster transitions. An example for wild-type ChR2 is given in Figure Figure4,4, where*t*_{IPI}100 ms.**Flux dependence**:*Off-curve*: {*G*_{d(1, 2)}, [*G*_{f0},*G*_{b0}]};*On-curve*: {All other parameters}. Voltage clamp (preferably): −70 mV, long pulse to steady-state, (e.g.,*T*≈ 500 ms) plus decay of off-curve. Vary light intensity from near threshold to saturation (e.g., ϕ = {0.1, 0.5, 1, 5, 10, 50, 100} mW/mm^{2},*n*≥ 5). The recorded off- and on-curves are automatically fitted. An example set is shown in Figure Figure55 with more details of the algorithm given in Appendix Section (Model-Dependent Fitting Procedures).

Additionally the six-state model requires one or more very short pulses in order to characterize the opsin activation rates which model the lag in transitioning to conductive states upon light stimulation:

**Opsin activation rate**: {*G*_{o1},*G*_{o2}} One or more*short*pulses, voltage clamp: −70 mV. Vary pulse length, e.g., 0.5, 1, 2, 3, varied up to 10 ms. PyRhO will automatically find the time of the peak current and use an iterative formula to estimate*G*_{o1}. We initially assume*G*_{o2}=*G*_{o1}. Further details of the algorithm are given in Appendix section Six-state opsin activation rate fitting (Step 1b.).

All light pulses should be “rectangular” (step functions) in that they have a sharp onset and offset. Examples of each protocol are included in PyRhO with illustrations provided in Figures Figures33–**7**. The duration of the on- and off-phases should also be kept approximately equal since the optimizer will effectively weight the contributions of each according to the relative numbers of data points. Additional parameter dependencies will be added in the future which may require additional data sets for a full characterization of the models.

In general, the most important data are those described for characterizing the flux dependence, which may be considered to be the “minimal set.” If this set consists of only a single photocurrent, the fitting algorithms will fix the parameters which model the flux dependence (ϕ_{m}, *p* and *q*) to the initial values supplied (along with fixing those describing other uncharacterized variables) and tune the remaining parameters to return a model fit for that specific flux. This is not recommended however, as the model is under-constrained by the data (typically resulting in a poorer fit than when using a whole set of photocurrents) and is unlikely to generalize well to new experimental conditions. For best results, the flux dependence photocurrents should be measured at light intensities spanning several orders of magnitude as described above.

If variations in other parameters or short pulses are of interest then the additional data should (ideally) be collected as described. However, if obtaining the data for a full characterization of the model is not possible, the pre-set default values should be adequate for most practical purposes.

Each voltage-clamp recorded photocurrent should be loaded into a `PhotoCurrent` object as follows:

pc = PhotoCurrent(I=i0, t=t, pulses=[[t_on, t_off]], phi=2e15, V=-70)

Here, `I` is the array of photocurrent values in nanoamperes, `t` is the corresponding array of recording times (or a scalar representing the time-step) in milliseconds, `pulses` is a nested list of *n* lists (or an *n* × 2 array), where *n* corresponds to the number of pulses and each inner list contains the on-time and off-time for each pulse in milliseconds, `phi` represents the stimulating flux value in photons · mm^{−2} · s^{−1} and `V` is the clamp voltage in millivolts (or “`None`” if the voltage was not clamped).

The `PhotoCurrent` class contains methods which automatically check the data and extract the key features from it, which may then be accessed as properties of the object with the `.` operator. Any properties which are data-derived are suffixed with “`_`,” for example, the peak and steady-state current are accessed with `pc.peak_` and `pc.ss_`, respectively. These photocurrents may easily be plotted, along with their main features and the light stimulus using the `plot()` method.

The `PhotoCurrent` objects are then combined into `ProtocolData` objects. These sets of photocurrents also provide several convenient methods for plotting and extracting parameters from the data set as a whole.

stepPD = ProtocolData(protocol="step", nRuns=1, phis=[1e14,1e15,1e16,1e17,1e18], Vs=[-70]) for iPhi, phi in enumerate(phis): for iV, V in enumerate(Vs): pc = PhotoCurrent(Is[iPhi][iV], t, pulses, phi, V) stepPD.addTrial(pc)

Finally, the data sets are combined into a dictionary using the protocol names as keys:

ChR2DataSet = {"step" : stepPD, "recovery" : recovPD, "rectifier" : rectiPD, "shortPulse" : shortPD}

This dictionary contains all the data necessary to parameterize all three models, however, if only the three and four-state models are of interest then the `"shortPulse"` protocol may be omitted.

Once the data have been loaded into the appropriate structures, the fitting algorithms may be called with the `fitModels()` function.

fp = fitModels(ChR2dataSet, nStates=6, params=initialParams)

This procedure returns a `Parameters` object (from the `lmfit` module) with the calculated values and plots the resultant model fits over the experimental photocurrents. The entire set of ChR2 data are shown fitted to the six-state model for each flux value spanning two orders of magnitude (ϕ = [2.21 × 10^{15}, 2.65 × 10^{17}] photons · mm^{−2} · s^{−1}) with the same set of parameters in Figure Figure5.5. The lowest and highest intensity photocurrents are also shown in Figure Figure66 with the model fits for both the three- and four-state models for direct comparison. The six-state model fits are not replotted here as they only exhibit a significant difference to the four-state fits for short pulses, as illustrated in Figure Figure77.

Having fit a model, it may be easily characterized by plotting how the light-dependent transition rates vary as a function of flux (based on the Hill equation) along with light-independent transition rates as shown in Figure Figure8.8. An individual fit is shown in more detail with the residual error for ϕ = 2.21 × 10^{15} photons · mm^{−2} · s^{−1} in Figure Figure9.9. To provide insight into the model's kinetics, PyRhO also offers state variable plots. The evolution of the six-state model corresponding to the fit in Figure Figure99 is given in Figure Figure1010.

In general terms, the fitting algorithm first finds the model-independent variables such as the dark recovery rate and voltage dependence factors, proceeding through “off-curve” parameters by fitting a double exponential decay function, optionally fitting opsin activation rates for the six-state model and finally optimizing across a set of “on-curves” to find any remaining parameters. Due to the inherent variability and imprecision in experimental measurements there is an optional second optimization phase over the entire set of photocurrents simultaneously. The values found for the dark parameters {*G*_{d(1, 2)}, [*G*_{f0}, *G*_{b0}]} (and opsin activation rates *G*_{o(1, 2)} if relevant) are used as the initial values, lower and upper bounds are calculated as 50 and 200% of these values, respectively (set by a hyperparameter) and the model is then re-optimized to achieve an overall better fit. The main sub-routines of the algorithm are given in Algorithm Algorithm11 with more detail for each process given in the Appendix.

When fitting the three-state model, a double exponential is fit (with two corresponding decay rates *G*_{d1} and *G*_{d2}) which are then weighted by their coefficients (*I*_{slow} and *I*_{fast}) and combined to form a single exponential. The mean of these values is then calculated across a set of *N* photocurrents (Equation 7) and this value is then used in subsequent parts of the fitting algorithm.

$$\begin{array}{r}{G}_{d}=\frac{1}{N}{\displaystyle \sum _{n=1}^{N}}\frac{{I}_{{slow}_{n}}\times {G}_{d{1}_{n}}+{I}_{{fast}_{n}}\times {G}_{d{2}_{n}}}{{I}_{{slow}_{n}}+{I}_{{fast}_{n}}}\text{}.\end{array}$$

(7)

In order to test the algorithms, synthetic data was generated using the parameter values derived from fitting the six-state model to the ChR2 experimental data. The fitting procedure was then applied to these synthetic photocurrents to compare the newly derived parameters with the known values used to generate the synthetic data. The results of these two fitting processes using the “`powell`” optimization algorithm are shown in Table Table33 (along with the values used as initial estimates).

We first note that many of the computed parameters are very close to the true (original) values (all but two are within ±5% of the original values), especially in the context of the degree of experimental noise and measurement error which would typically accompany recordings of real neuronal data. One notable exception is *G*_{o2} which is hard to fit in a single-pulse protocol since all opsin models are assumed to be in their fully dark-adapted (ground) state i.e., *C*_{1} = 1.

While there are other differences between some of the original and computed parameters, these may potentially be accounted for by numerical precision issues and the high-dimensional parameter space of the six-state model being under-constrained by the data. For example, a decrease in one parameter may be compensated for by an increase in another such that only latent variables are affected and the fit in the observable current is still good. Inspecting the model fit plots appears to confirm this, as the residual error is very low across the whole set of generated verification photocurrents—at most ±0.5% of the steady-state current and usually considerably less. The entire set of photocurrents fitted to the synthetic data are plotted in Figure Figure1111 and show a very close correspondence to the target (synthetic) photocurrents across the whole set of stimulus intensities.

To programmatically simulate the opsin, a model, protocol and simulator object must first be created (see Figure Figure1).1). The model and protocol are then loaded into the simulator which configures the simulation environment for that particular choice of opsin and protocol (e.g., by setting the numerical time-step according to the shortest stimulation period). The simulator object can then be run and plotted with the appropriate methods as shown in the following example, where a six-state model is used with the ChR2 parameters and run through the rectifier protocol (described below). In this example the protocol is run on the NEURON simulator, (rather than the default Python simulator) and the default ChR2 parameters are loaded from the module's `modelFits` dictionary, which contains some pre-fit model parameter sets for several common opsins.

from pyrho import * nStates = "6" # 3, 4 or 6 ChR2params = modelFits[nStates]["ChR2"] RhO = models[nStates](ChR2params) Prot = protocols["rectifier"]() Sim = simulators["NEURON"](Prot, RhO) Sim.run() Sim.plot()

Alternatively, when using the PyRhO GUI, the model, protocol and simulator are simply selected from drop-down lists, (optionally parameters may be changed), then simulated and plotted by clicking the “Run” button.

Each type of object will be initialized with default parameters (in the form of `Parameters` objects) unless passed a different set upon initialization e.g., `RhO =`
`models[nStates](params6s)`. Alternatively, parameters may be set after creation using methods such as `.setParams()` or `.updateParams()` for partial sets.

PyRhO comes with several preconfigured and customizable simulation protocols for exploring the dynamics of the models. These include typical system analysis stimuli such as delta functions, step functions and sinusoids, as well as chirps and specialized protocols designed to probe particular features of the opsins including voltage-dependence (`rectifier`), opsin activation (`shortPulse`) and dark recovery (`recovery`).

PyRhO's simulation layer serves to perform house-keeping tasks necessary to prepare different simulation environments to use a particular opsin model in a “system” of interest and apply a particular protocol to it. Currently, three simulators are available in PyRhO: Brian2 for neural networks, NEURON for detailed morphological neurons and (pure) Python for basic opsin channel dynamics. This selection of simulators provides PyRhO with a way to seamlessly span multiple scales of modeling with the same parameterized opsins; from individual channels to whole brain regions.

When using simulators other than “`Python`,” additional parameters may be specified. For example, the `NEURON` simulator has additional parameters such as “`v_init`” (the cell membrane potential initialization value), “`CVode`” (a boolean value for activating variable time-step solvers), and “`cell`” (a hoc file specifying the neuron's morphology). This allows existing simulations created in `NEURON` to be conveniently transfected (augmented with opsins) and run within the PyRhO framework. While models may be fit to data and then seamlessly inserted into one of these simulators within the same environment, if desired the `NMODL` files and Brian equations may be accessed and exported as a starting point for creating stand-alone simulations.

An example of implementing opsins within the `NEURON` environment is shown in Figure Figure1212 where ChR2 expressing cells are illuminated with a 150 ms light pulse. The six-state equations were used to model the ChR2 additions and implemented via the `NMODL` file `RhO6.mod`, adapted from the description in Grossman et al. (2013).

PyRhO also incorporates the Brian2 spiking neural network simulator. The opsin is represented with a set of ODEs which use the parameters specified in the `RhodopsinModel` object. As an example we show simulation results for a neural network which consists of 140 leaky integrate-and-fire neurons, separated into three feed-forward layers. The first group has neurons which express ChR2 and that layer has a set of random connections with the next layer, with 20% connectivity and 1*ms* conduction delays (with the same specifications for connectivity between the second and third layers). Figure Figure1313 shows raster plots of the spiking neurons in all three layers.

To make PyRhO usable without any programming background, a graphical user interface (GUI) has been written which runs in Jupyter (formerly IPython; Pérez and Granger, 2007). This runs in a browser-based notebook meaning that it could be easily configured as a server and made accessible to an entire laboratory or classroom without requiring local installations on each machine. Since both figures and text results appear embedded in the notebook after the code used to produce them, this makes the interface self-documenting and a particularly useful means of sharing models (Topalidou et al., 2015).

A screenshot of the GUI is given in Figure Figure1414 showing the simulation tab for the three-state model. In addition to the parameter fields, the model states diagram is also embedded in the GUI with the equations which describe the opsin's behavior rendered in LaTeX. Similarly, the parameters for each protocol and each simulator are shown on other tabs for easy modification of their values. On the fitting tab, there are also sub-tabs for each of the models, with their respective sets of parameters including tick-boxes to fix parameters and fields for numerical bounds and algebraic constraints to be passed to the fitting algorithms.

PyRhO has been written to be an integrated suite of intuitive, flexible, open-source and multi-scale computational tools for analysing and simulating opsins. In keeping with the open-source community's ethos, it builds upon existing libraries for numerical (NumPy), scientific (SciPy, lmfit) and plotting routines (Matplotlib). It also incorporates several freely available simulators, namely NEURON and Brian2, with work to incorporate PyNEST underway.

While the models were developed with rhodopsins and the fitting and simulation demonstrations are illustrated using ChR2 data, in principle the tools should work just as well with other classes of opsin. Essentially the opsin models represent non-linear dynamical systems (second order for the three-state model and *n* − 1^{th} order for the *n*-state model in general). As such they are capable of capturing the three main classes of dynamics in response to a step input: under-damped, over-damped and oscillatory (ringing) currents. While PyRhO can fit and simulate all three cases, interestingly, only the first two types of response have been observed so far (for low and high intensity stimulation, respectively).

For inhibitory opsins observed so far, while the sign of the current changes, the dynamics remain qualitatively the same, meaning that the fitting and simulation routines should be just as effective. However, for fitting, it may be necessary to adjust the starting values for some parameters to help the algorithms achieve good results. The main changes we anticipate would be in parameters that are extraneous to the the core dynamics (that is, those not described by the differential equations) such as the reversal potential, (*E*) which may be measured through standard electro-physiological techniques and possibly the parameters tuning the voltage rectification curve (*v*_{0} and *v*_{1}).

Given the way that PyRhO has been constructed around several layers of abstraction, extending the capabilities with more models, simulators or protocols is relatively straightforward. Nested classes (and sub-classes) act as templates guiding development while inheritance of methods and attributes facilitates code reuse in line with open-source software development best practices. As more data is collected, the library of parameterized models may be contributed to, making more characterized opsin variants available for other optogeneticists to simulate.

In broader terms, PyRhO has potential as a framework to be generalized to incorporate other types of stimulation, including for example electrical and magnetic. The simulator and protocol layers could potentially remain largely unaltered allowing it to relatively easily grow into a neural stimulation platform with neuro-engineering applications beyond the scope of just optogenetics.

The fitting algorithms rely upon optimization methods to extract several parameters and hence are subject to the standard issues associated with such procedures, including sensitivity to initial conditions and settling in local minima. To ameliorate these issues, the `lmfit` module is used to provide several useful facilities including imposing bounds on the parameters and algebraic constraints. Individual parameters may also be fixed and the fitting procedures rerun to optimize over the remaining free parameters, possibly with a new set of initial conditions. Measuring and then fixing physiological parameters in this way, such as the reversal potential *E*, will help the optimization algorithms and considerably improve the resultant model fit.

Additionally there are limitations due to the more specific nature of the models used. For example, the routine for estimating *g*_{0}, the biological scaling parameter for the cell's conductance, is systematically under-estimated. This should be the maximum conductance of the cell (voltage clamped to −70 mV) assuming full occupancy of the (primary) open-state. However, in simulations this condition is never achieved in the models, with maximum occupancy dependent on the model parameters, typically found to be around 0.8 for wild-type ChR2 (Nikolic et al., 2009). We therefore calculate a conservative correction factor to yield a better (albeit imperfect) estimate and make the user aware of the issue so that they may override and fix the values returned from the fitting function as necessary.

We note however that these limitations are not major, especially in the context of experimental and measurement inaccuracies and so do not significantly detract from the usefulness of PyRhO's fitting algorithms.

One natural development for PyRhO would be to extend the models (and fitting algorithms) to include additional parameter dependencies defined in the introduction, Equation (3), such as spectral absorption *f*_{λ}(λ), temperature *f*_{T}(*T*|*Q*_{10}), and pH *f*_{pH}(*pH*_{int}, *pH*_{ext}) factors for the channel conductance, as well as the effects of temperature and pH on the photocycle kinetics. This would allow the simulations to capture more of the variation and enhance the tools' ability to engineer the response to a target outcome or compensate for adverse experimental conditions.

The unconstrained nature of the parameter *G*_{o2} also suggests that there may be scope for additional models. For example, a simpler five-state model (similar to the six-state model but without *I*_{2} and its associated transitions) may be able to capture the experimental data as well as the six-state model with a less arduous fitting process and more stable resulting parameters (less sensitive to initial choices). In the first release however, we have used the six-state model for the sake of greater generality.

Another route for future development would be to incorporate other simulators for example PyNEST (Eppler et al., 2008), allowing greater flexibility for users to choose the simulator they are most comfortable with or which offers alternative features and performance benefits. Ultimately this process could be continued to interface with PyNN (Davison et al., 2009) and other simulator-agnostic model description languages such as NeuroML (Gleeson et al., 2010) and NineML (Gorchetchnikov et al., 2011). Furthermore, PyRhO could be combined with the software for optical pattern generation and data acquisition NeuroPG (Avants et al., 2015), to create a complete neural engineering tool for optogenetics.

While there is always scope to add additional protocols, a particularly interesting approach may be to allow networks of neurons to be transfected with several types of opsin and stimulated with multiple wavelengths of light (Han and Boyden, 2007). This could also be supplemented with a more detailed model of the optics of light-tissue interactions as previously implemented for an individual cell in NEURON (Foutz et al., 2012) or interfaced with OptogenSIM for whole brain simulations (Liu et al., 2015). These would be particularly useful additions to PyRhO's neuro-engineering capabilities, allowing stimuli to be more accurately sculpted. In the meantime, users of PyRhO should be mindful of the attenuating effects of light scattering and absorption (or equivalently a spectral shift) from passing through other tissue (particularly *in vivo*) which are not currently explicitly accounted for. These effects may result in a lower effective flux intensity than specified which may shift the “true” value of ϕ_{m} recovered from the fitting procedures.

In summary, we have presented and verified a new integrated suite of open-source computational tools for optogenetics. PyRhO has been demonstrated to characterize opsins from experimental photocurrents, fitting kinetic model parameters to yield a functional understanding, helping to guide opsin choice and development. These models have been demonstrated in simulations across multiple scales, from channels to networks, by harnessing popular simulators such as NEURON and Brian2. PyRhO is also provided with a Jupyter browser-based GUI to facilitate its use and aid in model sharing. We have outlined some of its chief strengths along with its limitations and plans for future improvements. By releasing these tools as open-source, we hope that other computational neuroscientists will contribute features and expertise, accelerating progress in the rapidly growing field of optogenetics.

BE designed and wrote the software module and GUI. KN and BE developed the opsin models. BE and KN developed the model fitting algorithms. BE primarily wrote the article, ran the tests, analyzed the data and plotted the figures with participation from KN. BE, KN, SJ, and SS contributed to the project concepts, manuscript comments and discussion of results.

This work was supported by the UK Biotechnology and Biological Sciences Research Council (BBSRC) grant BB/L018268/1 and the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/N002474/1.

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.

We would like to thank Matthew Grub and Juan Burrone for the ChR2 photocurrent data used to illustrate the fitting algorithms. We would also like to thank Pepe Herrero for designing the PyRhO logo.

The analytic solution for arbitrary initial conditions (*C*_{0}, *O*_{0}, *D*_{0}) is as follows:

$$\begin{array}{l}C(t)=\frac{{G}_{d}{G}_{r}}{{\lambda}_{1}{\lambda}_{2}}+\frac{({\lambda}_{1}-{G}_{d}-{G}_{r})}{{G}_{d}{\lambda}_{1}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{1}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{1}\xb7t}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}-\text{\hspace{0.17em}}\frac{({\lambda}_{2}-{G}_{d}-{G}_{r})}{{G}_{d}{\lambda}_{2}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{2}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{2}\xb7t}\end{array}$$

(A1)

$$\begin{array}{l}O(t)=\frac{{G}_{a}{G}_{r}}{{\lambda}_{1}{\lambda}_{2}}-\frac{({\lambda}_{1}-{G}_{r})}{{G}_{d}{\lambda}_{1}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{1}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{1}\xb7t}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}-\text{\hspace{0.17em}}\frac{({\lambda}_{2}-{G}_{r})}{{G}_{d}{\lambda}_{2}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{2}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{2}\xb7t}\end{array}$$

(A2)

$$\begin{array}{l}D(t)=\frac{{G}_{a}{G}_{d}}{{\lambda}_{1}{\lambda}_{2}}+\frac{1}{{\lambda}_{1}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{1}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{1}\xb7t}\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}-\text{\hspace{0.17em}}\frac{1}{{\lambda}_{2}\xi}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{\zeta}_{2}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{e}^{-{\lambda}_{2}\xb7t}\end{array}$$

(A3)

with the following definitions:

$$\xi =\sqrt{{G}_{a}^{2}+{G}_{d}^{2}+{G}_{r}^{2}-2\xb7({G}_{a}{G}_{d}+{G}_{a}{G}_{r}+{G}_{d}{G}_{r})}$$

$${\lambda}_{1,2}=\frac{1}{2}\xb7\left(({G}_{a}+{G}_{d}+{G}_{r})\pm \xi \right)$$

$$\begin{array}{rcl}{\zeta}_{1,2}& =& {C}_{0}\xb7\left[{G}_{d}\xb7{G}_{a}\right]+{O}_{0}\xb7\left[{G}_{d}\xb7({G}_{a}-{\lambda}_{1,2})\right]\end{array}$$

$$\begin{array}{r}+\text{}{D}_{0}\xb7\left[{G}_{r}\xb7({\lambda}_{1,2}-{G}_{a}-{G}_{d})\right]\end{array}$$

Note also that λ_{1} · λ_{2} = *G*_{a}*G*_{d} + *G*_{a}*G*_{r} + *G*_{d}*G*_{r} i.e., the sum of products.

Prior to fitting the specific parameters of a chosen model, several fitting processes are performed to find parameters common to all models:

- If a set of voltage clamp data are provided over a range of potentials, the steady-state currents are used to fit the rectification function,
*f*_{v}(*v*), (Equation 5) yielding the parameters*E*,*v*_{0}, and*v*_{1}:- The
*I*_{ss}values are first fit to the function*f*_{v}(**V**) · (_{clamp}**V**−_{clamp}*E*) to find an initial estimate for*v*_{0}and a value for*E*where the polarity of the steady-state currents changes (along with a dummy parameter comprising*g*_{0}·*f*_{ϕ}·*v*_{1}, which is discarded). - The reversal potential,
*E*, is fixed and the steady-state currents are converted to conductances:**g**=**I**−_{ss}/ (V_{clamp}*E*). These conductances are then normalized by dividing by the conductance calculated at*V*_{clamp}= −70 mV to produce an array of conductance factors. - The conductance factors are fit to the function
*f*_{v}(*V*_{clamp}) with the constraint ${v}_{1}=\left(\mathrm{70}+E\right)\u2215\left({e}^{\frac{70+E}{{v}_{0}}}-1\right)$ to calculate and fix the parameters*E*,*v*_{0}and*v*_{1}.

- All photocurrents in the dataset recorded at −70
*mV*are scanned to find the highest peak. This value is used to calculate an initial estimate for*g*_{0}(according to the formula*g*_{0}=*I*_{peak, max}/ (−70 −*E*)) to aid the optimization routines. - If the recovery protocol is included in the dataset (pairs of pulses with a range of inter-pulse-intervals) then the peak currents of the second pulses are extracted along with the times and the initial peak and steady-state current. These values are then aligned to the end of the first pulse (where the current is
*I*_{ss0}) and fit to the exponential recovery function: ${I}_{\text{peak}}(t)={I}_{peak0}-a\xb7{e}^{{G}_{r0}\xb7t}$. This yields the dark recovery rate constant,*G*_{r0}.

These parameters are then fixed (except for the estimate of *g*_{0}) and passed to the model-specific fitting routines.

Each model is fit with broadly the same three-stage fitting procedure consisting of Off-curve parameters, then On-curve parameters and finally re-optimizing:

- 1a. Fit bi-exponential functions to the Off-curves to find parameter values for the light-insensitive decay rates: {
*G*_{d1},*G*_{d2},*G*_{f0},*G*_{b0}} for the four- and six-state models [detailed in Section Four- and Six-state Model Off-curve Fitting (Step 1a.)], or*G*_{d}for the three-state model (described in Equation 7), - 1b. [Six-state only] Fix the values for the parameters from Step 1a then calculate activation rates
*G*_{o1}and*G*_{o2}[detailed in Section Six-State Opsin Activation Rate Fitting (Step 1b.)], - 2. Fix the values for the parameters found in Step 1 and fit the On-curves to find the remaining parameters governing the light-dependent transitions,
- 3. Using the values found in Steps 1 and 2 as initial values, relax all parameter values to freely vary within a range (by default between 50 and 200% of the initial values) and refit the model across the set of whole photocurrent curves (on- and off-phases).

The fitting equations for the Light-Off response curves can be calculated directly from the set of ODEs given in Table Table2,2, for light flux zero, as shown in Nikolic et al. (2009):

$$\begin{array}{rcl}{O}_{1\text{off}}& =& {A}_{11}exp(-{\Lambda}_{1}t)+{A}_{12}exp(-{\Lambda}_{2}t)\end{array}$$

(A4)

$$\begin{array}{rcl}{O}_{2\text{off}}& =& {A}_{21}exp(-{\Lambda}_{1}t)+{A}_{22}exp(-{\Lambda}_{2}t)\end{array}$$

(A5)

where

$${\Lambda}_{1,2}=b\mp c\text{},$$

*b* = (*G*_{d1} + *G*_{d2} + *G*_{f0} + *G*_{b0})/2 and $c=\sqrt{{b}^{2}-({G}_{\text{d}1}{G}_{\text{d}2}+{G}_{\text{d}1}{G}_{\text{b}0}+{G}_{\text{d}2}{G}_{\text{f}0})}$. Two useful relationships follow from the previous equations:

$$\begin{array}{rcl}{\Lambda}_{1}+{\Lambda}_{2}& =& {G}_{\text{d}1}+{G}_{\text{d}2}+{G}_{\text{f}0}+{G}_{\text{b}0}\end{array}$$

(A6)

$$\begin{array}{rcl}{\Lambda}_{1}\xb7{\Lambda}_{2}& =& {G}_{\text{d}1}{G}_{\text{d}2}+{G}_{\text{d}1}{G}_{\text{b}0}+{G}_{\text{d}2}{G}_{\text{f}0}\end{array}$$

(A7)

The current after the light is turned off is given by:

$$\begin{array}{l}{I}_{\text{off}}=V({g}_{1}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{O}_{1\text{off}}+{g}_{2}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{O}_{2\text{off}})\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}={I}_{\text{slow}}\mathrm{exp}(-{\Lambda}_{1}t)\text{\hspace{0.17em}}+{I}_{\text{fast}}\mathrm{exp}(-{\Lambda}_{2}t),\end{array}$$

(A8)

where the expression for *I*_{slow} (corresponds to Λ_{1}) and *I*_{fast} (corresponds to Λ_{2}) components of the off-current can be calculated from the equation above and are given in Nikolic et al. (2009). Note also that *I*_{slow} + *I*_{fast} = *I*_{ss} which is included as an additional constraint in the fitting procedure.

From the differential equations describing the dynamics of a short pulse (Table (Table2,2, three-state model):

$$\frac{dO}{dt}={G}_{a}(t)C-{G}_{d}O$$

where *G*_{a}(*t*) is given in Nikolic et al., 2007, Appendix 1. If the light illumination occurs from *t* = 0 to *t* = *t*_{pulse}, then after the light goes off (*t* > *t*_{pulse}):

$$\begin{array}{l}{G}_{a}(t)={G}_{a}\left[{e}^{{t}_{pulse}\u2215{\tau}_{opsin}}-1\right]\xb7{e}^{-t\u2215{\tau}_{opsin}}\text{},\end{array}$$

(A9)

where τ_{opsin}(= 1/*G*_{o}) is the time constant of the opsin complex activation (or conformal change) after the retinal isomerization (which happens very quickly upon photon absorption). Now if we use a very short pulse then *t*_{pulse} τ_{opsin} so that *C* ≈ 1 (and secondary cycles are not significantly activated), we may apply the standard approximation for small exponents *e*^{x} ≈ 1 + *x*:

$${G}_{a}(t)\approx {G}_{a}\frac{{t}_{pulse}}{{\tau}_{opsin}}\xb7{e}^{-t\u2215{\tau}_{opsin}}=Q\xb7{e}^{-\frac{t}{{\tau}_{opsin}}}\text{},\text{where}Q={G}_{a}\frac{{t}_{pulse}}{{\tau}_{opsin}}$$

and so

$$\begin{array}{l}\frac{dO}{dt}\approx Q\xb7{e}^{-{G}_{o}t}-{G}_{d}O\text{}.\end{array}$$

(A10)

The solution of Equation (A10) is the limit case for very short pulses (relative to τ_{opsin}) when *I*(*t*) stops being dependent on the pulse duration *t*_{pulse}:

$$\begin{array}{l}I(t)=g\xb7(v-E)\xb7\frac{Q}{({G}_{o}-{G}_{d})}\xb7{e}^{-t\u2215{\tau}_{d}}\xb7(1-{e}^{-({G}_{o}-{G}_{d})t})\text{}.\end{array}$$

(A11)

Differentiating Equation (A11) we find the time of the peak current:

$${t}_{maxlag}=\frac{log({G}_{o}\u2215{G}_{d})}{{G}_{o}-{G}_{d}}=\frac{{\tau}_{opsin}{\tau}_{d}}{{\tau}_{d}-{\tau}_{opsin}}log\left(\frac{{\tau}_{d}}{{\tau}_{opsin}}\right)\text{},$$

where τ_{d} = 1/*G*_{d}. Assuming that τ_{opsin} τ_{d} then,

$$\begin{array}{l}{t}_{maxlag}\approx {\tau}_{opsin}log\left(\frac{{\tau}_{d}}{{\tau}_{opsin}}\right)\text{}.\end{array}$$

(A12)

To find *t*_{maxlag} we can fit the following function to a series of short pulses:

$${t}_{peak}={t}_{pulse}+{t}_{maxlag}\xb7{e}^{-k\xb7{t}_{pulse}}\text{}.$$

Alternatively, estimate *t*_{maxlag} from the “delta” protocol as *t*_{pulse} → 0. Once we have obtained an estimate for *t*_{maxlag} we can then solve for τ_{opsin} iteratively as follows:

$$\begin{array}{l}{\tau}_{opsin,i+1}=\frac{{t}_{maxlag}}{({t}_{maxlag}\u2215{\tau}_{d})-log({\tau}_{opsin,i}\u2215{\tau}_{d})}\text{}.\end{array}$$

(A13)

- Adamantidis A. R., Zhang F., Aravanis A. M., Deisseroth K., de Lecea L. (2007). Neural substrates of awakening probed with optogenetic control of hypocretin neurons. Nature 450, 420–424. 10.1038/nature06310 [PubMed] [Cross Ref]
- Airan R. D., Thompson K. R., Fenno L. E., Bernstein H., Deisseroth K. (2009). Temporally precise
*in vivo*control of intracellular signalling. Nature 458, 1025–1029. 10.1038/nature07926 [PubMed] [Cross Ref] - 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. 10.1088/1741-2560/4/3/s02 [PubMed] [Cross Ref] - Arenkiel B. R., Peca J., Davison I. G., Feliciano C., Deisseroth K., Augustine G. J., et al. . (2007).
*In vivo*light-induced activation of neural circuitry in transgenic mice expressing channelrhodopsin-2. Neuron 54, 205–218. 10.1016/j.neuron.2007.03.005 [PMC free article] [PubMed] [Cross Ref] - Arlow R. L., Foutz T. J., McIntyre C. C. (2013). Theoretical principles underlying optical stimulation of myelinated axons expressing channelrhodopsin-2. Neuroscience 248, 541–551. 10.1016/j.neuroscience.2013.06.031 [PMC free article] [PubMed] [Cross Ref]
- Arrenberg A. B., Stainier D. Y. R., Baier H., Huisken J. (2010). Optogenetic control of cardiac function. Science 330, 971–974. 10.1126/science.1195929 [PubMed] [Cross Ref]
- Avants B. W., Murphy D. B., Dapello J. A., Robinson J. T. (2015). Neuropg: open source software for optical pattern generation and data acquisition. Front. Neuroeng. 8:1 10.3389/fneng.2015.00001 [PMC free article] [PubMed] [Cross Ref]
- Ayling O. G. S., Harrison T. C., Boyd J. D., Goroshkov A., Murphy T. H. (2009). Automated light-based mapping of motor cortex by photoactivation of channelrhodopsin-2 transgenic mice. Nat. Methods 6, 219–224. 10.1038/nmeth.1303 [PubMed] [Cross Ref]
- 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. 10.1038/ncomms6810 [PubMed] [Cross Ref]
- Bamann C., Kirsch T., Nagel G., Bamberg E. (2008). Spectral characteristics of the photocycle of channelrhodopsin-2 and its implication for channel function. J. Mol. Biol. 375, 686–694. 10.1016/j.jmb.2007.10.072 [PubMed] [Cross Ref]
- Bamberg E., Bamann C., Feldbauer K., Kleinlogel S., Spitz J., Zimmermann D., et al. (2008). Channelrhodopsins: Molecular Properties and Applications. Washington, DC: Society for Neuroscience.
- Berndt A., Yizhar O., Gunaydin L. A., Hegemann P., Deisseroth K. (2008). Bi-stable neural state switches. Nat. Neurosci. 12, 229–234. 10.1038/nn.2247 [PubMed] [Cross Ref]
- Boyden E. S., Zhang F., Bamberg E., Nagel G., Deisseroth K. (2005). Millisecond-timescale, genetically targeted optical control of neural activity. Nat. Neurosci. 8, 1263–1268. 10.1038/nn1525 [PubMed] [Cross Ref]
- Boyle P. M., Williams J. C., Ambrosi C. M., Entcheva E., Trayanova N. A. (2013). A comprehensive multiscale framework for simulating optogenetics in the heart. Nat. Commun. 4, 1–9. 10.1038/ncomms3370 [PMC free article] [PubMed] [Cross Ref]
- Bruegmann T., Sasse P. (2015). Optogenetic cardiac pacemakers: science or fiction? Trends Cardiovas. Med. 25, 82–83. 10.1016/j.tcm.2014.10.016 [PubMed] [Cross Ref]
- Busskamp V., Roska B. (2011). Optogenetic approaches to restoring visual function in retinitis pigmentosa. Curr. Opin. Neurobiol. 21, 1–5. 10.1016/j.conb.2011.06.001 [PubMed] [Cross Ref]
- Chow B. Y., Han X., Dobry A. S., Qian X., Chuong A. S., Li M., et al. . (2010). High-performance genetically targetable optical neural silencing by light-driven proton pumps. Nature 463, 98–102. 10.1038/nature08652 [PMC free article] [PubMed] [Cross Ref]
- Chuong A. S., Miri M. L., Busskamp V., Matthews G. A. C., Acker L. C., Sørensen A. T., et al. . (2014). Noninvasive optical inhibition with a red-shifted microbial rhodopsin. Nat. Neurosci. 17, 1123–1129. 10.1038/nn.3752 [PMC free article] [PubMed] [Cross Ref]
- Davison A. P., Brüderle D., Eppler J. M., Kremkow J., Muller E., Pecevski D., et al. (2009). Pynn: a common interface for neuronal network simulators. Front. Neuroinform. 2:11 10.3389/neuro.11.011.2008 [PMC free article] [PubMed] [Cross Ref]
- Degenaar P., Grossman N., Memon M. A., Burrone J., Dawson M., Drakakis E., et al. . (2009). Optobionic vision: a new genetically enhanced light on retinal prosthesis. J. Neural Eng. 6:035007. 10.1088/1741-2560/6/3/035007 [PubMed] [Cross Ref]
- Eppler J. M., Helias M., Muller E., Diesmann M., Gewaltig M.-O. (2008). Pynest: A convenient interface to the nest simulator. Front. Neuroinform. 2:12 10.3389/neuro.11.012.2008 [PMC free article] [PubMed] [Cross Ref]
- Ernst O. P., Sánchez Murcia P. A., Daldrop P., Tsunoda S. P., Kateriya S., Hegemann P. (2008). Photoactivation of channelrhodopsin. J. Biol. Chem. 283, 1637–1643. 10.1074/jbc.M708039200 [PubMed] [Cross Ref]
- Feldbauer K., Zimmermann D., Pintschovius V., Spitz J., Bamann C., Bamberg E. (2009). Channelrhodopsin-2 is a leaky proton pump. Proc. Natl. Acad. Sci. U.S.A. 106, 12317–12322. 10.1073/pnas.0905852106 [PubMed] [Cross Ref]
- Foutz T. J., Arlow R. L., McIntyre C. C. (2012). Theoretical principles underlying optical stimulation of a channelrhodopsin-2 positive pyramidal neuron. J. Neurophysiol. 107, 3235–3245. 10.1152/jn.00501.2011 [PubMed] [Cross Ref]
- Gleeson P., Crook S., Cannon R. C., Hines M. L., Billings G. O., Farinella M., et al. (2010). Neuroml: a language for describing data driven models of neurons and networks with a high degree of biological detail. PLoS Comput. Biol. 6:e1000815 10.1371/journal.pcbi.1000815 [PMC free article] [PubMed] [Cross Ref]
- Goodman D., Brette R. (2008). Brian: a simulator for spiking neural networks in python. Front. Neuroinform. 2:5. 10.3389/neuro.11.005.2008 [PMC free article] [PubMed] [Cross Ref]
- Goodman D. F. M., Brette R. (2009). The brian simulator. Front. Neurosci. 3, 192–197. 10.3389/neuro.01.026.2009 [PMC free article] [PubMed] [Cross Ref]
- Gorchetchnikov A., Cannon R., Clewley R., Cornelis H., Davison A., De Schutter E., et al. (2011). Nine
*ML*: declarative, mathematically-explicit descriptions of spiking neuronal networks. Front. Neuroinform. Conference Abstract: 4th INCF Congress of Neuroinformatics. 10.3389/conf.fninf.2011.08.00098 [Cross Ref] - Gradinaru V., Mogri M., Thompson K. R., Henderson J. M., Deisseroth K. (2009). Optical deconstruction of parkinsonian neural circuitry. Science 324, 354–359. 10.1126/science.1167093 [PubMed] [Cross Ref]
- Gradinaru V., Thompson K. R., Zhang F., Mogri M., Kay K., Schneider M. B., et al. . (2007). Targeting and readout strategies for fast optical neural control
*in vitro*and*in vivo*. J. Neurosci. 27, 14231–14238. 10.1523/JNEUROSCI.3578-07.2007 [PubMed] [Cross Ref] - Gradmann D., Berndt A., Schneider F., Hegemann P. (2011). Rectification of the channelrhodopsin early conductance. Biophys. J. 101, 1057–1068. 10.1016/j.bpj.2011.07.040 [PubMed] [Cross Ref]
- Gradmann D., Ehlenbeck S., Hegemann P. (2002). Modeling light-induced currents in the eye of chlamydomonas reinhardtii. J. Mem. Biol. 189, 93–104. 10.1007/s00232-002-1006-8 [PubMed] [Cross Ref]
- Grossman N., Nikolic K., Toumazou C., Degenaar P. (2011). Modeling study of the light stimulation of a neuron cell with channelrhodopsin-2 mutants. IEEE Trans. Biomed. Eng. 58, 1742–1751. 10.1109/TBME.2011.2114883 [PubMed] [Cross Ref]
- Grossman N., Simiaki V., Martinet C., Toumazou C., Schultz S. R., Nikolic K. (2013). The spatial pattern of light determines the kinetics and modulates backpropagation of optogenetic action potentials. J. Comp. Neurosci. 34, 477–488. 10.1007/s10827-012-0431-7 [PMC free article] [PubMed] [Cross Ref]
- Gunaydin L. A., Yizhar O., Berndt A., Sohal V. S., Deisseroth K., Hegemann P. (2010). Ultrafast optogenetic control. Nat. Neurosci. 13, 387–392. 10.1038/nn.2495 [PubMed] [Cross Ref]
- Han X., Boyden E. S. (2007). Multiple-color optical activation, silencing, and desynchronization of neural activity, with single-spike temporal resolution. PLoS ONE 2:e299. 10.1371/journal.pone.0000299 [PMC free article] [PubMed] [Cross Ref]
- Hegemann P., Ehlenbeck S., Gradmann D. (2005). Multiple photocycles of channelrhodopsin. Biophys. J. 89, 3911–3918. 10.1529/biophysj.105.069716 [PubMed] [Cross Ref]
- Hegemann P., Möglich A. (2011). Channelrhodopsin engineering and exploration of new optogenetic tools. Nat. Methods 8, 39–42. 10.1038/nmeth.f.327 [PubMed] [Cross Ref]
- Hernandez V. H., Gehrt A., Reuter K., Jing Z., Jeschke M., Mendoza Schulz A., et al. . (2014). Optogenetic stimulation of the auditory pathway. J. Clin. Invest. 124, 1114–1129. 10.1172/JCI69050 [PMC free article] [PubMed] [Cross Ref]
- Hines M., Davison A. P., Muller E. (2009). Neuron and python. Front. Neuroinform. 3:1. 10.3389/neuro.11.001.2009 [PMC free article] [PubMed] [Cross Ref]
- Hines M. L., Carnevale N. T. (2000). Expanding neuron's repertoire of mechanisms with nmodl. Neural Comput. 12, 995–1007. 10.1162/089976600300015475 [PubMed] [Cross Ref]
- Ishizuka T., Kakuda M., Araki R., Yawo H. (2006). Kinetic evaluation of photosensitivity in genetically engineered neurons expressing green algae light-gated channels. Neurosci. Res. 54, 85–94. 10.1016/j.neures.2005.10.009 [PubMed] [Cross Ref]
- 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. Meth. 11, 338–346. 10.1038/nmeth.2836 [PMC free article] [PubMed] [Cross Ref]
- Konermann S., Brigham M. D., Trevino A. E., Hsu P. D., Heidenreich M., Le C., et al. . (2013). Optical control of mammalian endogenous transcription and epigenetic states. Nature 500, 472–476. 10.1038/nature12466 [PMC free article] [PubMed] [Cross Ref]
- Kuhne J., Eisenhauer K., Ritter E., Hegemann P., Gerwert K., Bartl F. (2014). Early formation of the ion-conducting pore in channelrhodopsin-2. Angewandte Chemie 54, 4953–4957. 10.1002/anie.201410180 [PubMed] [Cross Ref]
- Lagali P. S., Balya D., Awatramani G. B., Münch T. A., Kim D. S., Busskamp V., et al. . (2008). Light-activated channels targeted to on bipolar cells restore visual function in retinal degeneration. Nat. Neurosci. 11, 667–675. 10.1038/nn.2117 [PubMed] [Cross Ref]
- Lin J. Y. (2011). A user's guide to channelrhodopsin variants: features, limitations and future developments. Exp. Physiol. 96, 19–25. 10.1113/expphysiol.2009.051961 [PMC free article] [PubMed] [Cross Ref]
- Lin J. Y., Lin M. Z., Steinbach P., Tsien R. Y. (2009). Characterization of engineered channelrhodopsin variants with improved properties and kinetics. Biophys. J. 96, 1803–1814. 10.1016/j.bpj.2008.11.034 [PubMed] [Cross Ref]
- Liu Y., Jacques S. L., Azimipour M., Rogers J. D., Pashaie R., Eliceiri K. W. (2015). Optogensim: a 3d monte carlo simulation platform for light delivery design in optogenetics. Biomed. Opt. Exp. 6, 4859–4870. 10.1364/BOE.6.004859 [PMC free article] [PubMed] [Cross Ref]
- Muller E., Bednar J. A., Diesmann M., Gewaltig M.-O., Hines M., Davison A. P. (2015). Python in neuroscience. Front. Neuroinform. 9:11. 10.3389/fninf.2015.00011 [PMC free article] [PubMed] [Cross Ref]
- 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. 10.1073/pnas.1936192100 [PubMed] [Cross Ref]
- Newville M., Stensitzki T., Allen D. B., Ingargiola A. (2014). LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python. Zenodo. 10.5281/zenodo.11813 [Cross Ref]
- Nikolic K., Grossman N., Grubb M. S., Burrone J., Toumazou C., Degenaar P. (2009). Photocycles of channelrhodopsin-2. Photochem. Photobiol. 85, 400–411. 10.1111/j.1751-1097.2008.00460.x [PubMed] [Cross Ref]
- Nikolic K., Grossman N., Yan H., Drakakis E., Toumazou C., Degenaar P. (2007). A non-invasive retinal prosthesis - testing the concept, in Engineering in Medicine and Biology Society, EMBS 2007. 29th Annual International Conference of the IEEE (Lyon: ), 6364–6367. [PubMed]
- Pérez F., Granger B. E. (2007). IPython: a system for interactive scientific computing. Comput. Sci. Eng. 9, 21–29. 10.1109/MCSE.2007.53 [Cross Ref]
- Petreanu L., Mao T., Sternson S. M., Svoboda K. (2009). The subcellular organization of neocortical excitatory connections. Nature 457, 1142–1145. 10.1038/nature07709 [PMC free article] [PubMed] [Cross Ref]
- Shoham S., Deisseroth K. (2010). Special issue on optical neural engineering: advances in optical stimulation technology. J. Neural Eng. 7:040201. 10.1088/1741-2560/7/4/040201 [PubMed] [Cross Ref]
- Stehfest K., Hegemann P. (2010). Evolution of the channelrhodopsin photocycle model. Chemphyschem 11, 1120–1126. 10.1002/cphc.200900980 [PubMed] [Cross Ref]
- Topalidou M., Leblois A., Boraud T., Rougier N. P. (2015). A long journey into reproducible computational neuroscience. Front. Comput. Neurosci. 9:30. 10.3389/fncom.2015.00030 [PMC free article] [PubMed] [Cross Ref]
- Wang H., Peca J., Matsuzaki M., Matsuzaki K., Noguchi J., Qiu L., et al. . (2007). High-speed mapping of synaptic connectivity using photostimulation in channelrhodopsin-2 transgenic mice. Proc. Natl. Acad. Sci. U.S.A. 104, 8143–8148. 10.1073/pnas.0700384104 [PubMed] [Cross Ref]
- 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. 10.1371/journal.pcbi.1003220 [PMC free article] [PubMed] [Cross Ref]
- Yizhar O., Fenno L.E., Davidson T. J, Mogri, M., Deisseroth K. (2011). Optogenetics in neural systems. Neuron 71, 9–34. 10.1016/j.neuron.2011.06.004 [PubMed] [Cross Ref]
- Zhang F., Vierock J., Yizhar O., Fenno L. E., Tsunoda S., Kianianmomeni A., et al. . (2011). The microbial opsin family of optogenetic tools. Cell 147, 1446–1457. 10.1016/j.cell.2011.12.004 [PMC free article] [PubMed] [Cross Ref]
- Zhang F., Wang L. P., Boyden E. S., Deisseroth K. (2006). Channelrhodopsin-2 and optical control of excitable cells. Nat. Methods 3, 785–792. 10.1038/nmeth936 [PubMed] [Cross Ref]
- Zhang F., Wang L. P., Brauner M., Liewald J. F., Kay K., Watzke N., et al. . (2007). Multimodal fast optical interrogation of neural circuitry. Nature 446, 633–639. 10.1038/nature05744 [PubMed] [Cross Ref]

Articles from Frontiers in Neuroinformatics are provided here courtesy of **Frontiers Media SA**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |