PyRhO: A Multiscale Optogenetics Simulation Platform

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.

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 opensource, 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: 1. The first tool will automatically fit a choice of models to experimental data, extracting the parameters that describe the functional dynamics of the opsins. 2. 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 voltageclamp 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. 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'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.

Implementation
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 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 1 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).

Photocurrent Model
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: 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): In this version of PyRhO we have implemented the photocycle and membrane voltage dependencies and assumed that these two contributions can be separated: 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.

Photocycle Models
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 1 with accompanying illustrations in Figure 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 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 2), with channel conductances g O 1 and g O 2 , respectively. The Photocycle factor in Equation (1) has the form: where O

Voltage Dependence
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 . 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: 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 Frontiers in Neuroinformatics | www.frontiersin.org which cause the exponent to be small and the expression in Equation (5) , 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: This relationship is used as a constraint in the fitting procedures described below.

Model Fitting
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.

Characterization Data
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 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 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 peak (t) = I peak0 −a·e −G r0 ·t 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 4, where t IPI 100 ms.
The recorded off-and on-curves are automatically fitted. An example set is shown in Figure 5 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 3-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.

Minimal Data Requirements
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.
FIGURE 5 | The six-state model fit to a set of six ChR2 photocurrents using the same model parameters.
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.
FIGURE 6 | The three-and four-state models fitted to the first two ChR2 photocurrents for comparison (the six-state fits are omitted since they are almost indistinguishable from the four-state fits at this time scale). While both models capture the major features of the data, the four-state (and six-state) models produce a higher quality of fit, particularly during the post-peak inactivation phase (where the three-state curves tend to over-shoot the data) and the deactivation (off) phases which are better described by a double rather than a single exponential decay process. 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.

Fitting Procedure and Algorithms
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 5. The lowest and highest intensity photocurrents are also shown in Figure 6 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 7.
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 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 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 9 is given in Figure 10.
In general terms, the fitting algorithm first finds the modelindependent 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 f 0 , 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 FIGURE 7 | The four-state and six-state models simulated with ChR2-derived parameters for short-duration pulses. With the four-state photocurrents (A), the peak occurs at the end of the on-phase. This can be seen most clearly in the plot of Pulse duration vs. Time of peak (B) where all points lie on the diagonal. For the six-state photocurrents however, the peaks lag behind the end of the illumination periods slightly (D), due to the transitions to and from the model's extra intermediate states. The inactivation phases (A,D) and magnitude of the peaks (C,F) can be seen to be very similar for both models. For the six-state model, the peak-lag effect can also be observed to diminish as the pulse duration increases (E), demonstrating the practical equivalence of the two models for long stimulation periods.
given in Algorithm 1 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.
I slow n × G d1 n + I fast n × G d2 n I slow n + I fast n .

Verification
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 3 (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 if "rectifier" in {DataSet} then 3: if "recovery" in {DataSet} then 6: G r0 ← fitPeakRecovery(t p , I p ) 7: (G d (1,2) , if postFitOptimization is True then 10: 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 11 and show a very close correspondence to the target (synthetic) photocurrents across the whole set of stimulus intensities.

Simulations Procedure
To programmatically simulate the opsin, a model, protocol and simulator object must first be created (see Figure 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. Sim.run() Sim.plot() Alternatively, when using the PyRhO GUI, the model, protocol and simulator are simply selected from dropdown lists, (optionally parameters may be changed), FIGURE 9 | An example six-state model fit to ChR2 data (at φ = 2.21 × 10 15 photons · mm −2 · s −1 ) with the accompanying residual error expressed as a percentage of the steady-state current.
FIGURE 10 | The six-state model's internal states. This figure corresponds to the example fitting plot (Figure 9) and shows the evolution of the internal states through time (A,B) shows alternative representations of the same data) along with the occupancy proportions at the initial conditions (C), peak (D) and steady-state (E).
Frontiers in Neuroinformatics | www.frontiersin.org 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.

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

Simulators
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 12 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 1ms conduction delays (with the same specifications for connectivity between the second and third layers). Figure 13 shows raster plots of the spiking neurons in all three layers.

Graphical User Interface
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 14 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.

DISCUSSION
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 opensource community's ethos, it builds upon existing libraries for numerical (NumPy), scientific (SciPy, lmfit) and plotting FIGURE 11 | The six-state model (re-)fitted to synthetic data generated from parameters derived from six ChR2 photocurrents. It can be seen that the model fits (black lines) almost perfectly fit each synthetic photocurrent using a unified parameter set. routines (Matplotlib). It also incorporates several freely available simulators, namely NEURON and Brian2, with work to incorporate PyNEST underway.

Classes of Opsin
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 nonlinear 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 ).

Extensibility
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.

Limitations
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.

Future Work
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 ) 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.

AUTHOR CONTRIBUTIONS
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.

FUNDING
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.