A Virtual Optogenetics Laboratory


In [ ]:
# Import the module and set figure rendering to be inline
%matplotlib inline
from pyrho import *
setFigOutput('screen')

Graphical User Interface


In [ ]:
loadGUI()

Example fitting process


In [ ]:
initParams = Parameters()
initParams.add_many(                
                # Name   Value   Vary    Min     Max     Expr
                ('g0',   2.5e4,  True,   0.0,    1e15,   None),
                ('gam',  0.05,   True,   0.00,   1,      None),
                ('phi_m',3.5e17, True,   1e15,   1e19,   None),
                ('k1',   10,     True,   0.0,    1000,   None),
                ('k2',   3,      True,   0.0,    1000,   None),
                ('p',    1,      True,   0.1,    5,      None),
                ('Gf0',  0.04,   True,   0.0,    1000,   None),
                ('k_f',  0.1,    True,   0.0,    1000,   None),
                ('Gb0',  0.02,   True,   0.0,    1000,   None),
                ('k_b',  0.15,   True,   0.0,    1000,   None),
                ('q',    1,      True,   0.1,    5,      None),
                ('Go1',  2,      True,   0.0,    1000,   None),
                ('Go2',  2,      True,   0.0,    1000,   None),
                ('Gd1',  0.1,    True,   0.0,    1000,   None),
                ('Gd2',  0.01,   True,   0.0,    1000,   None),
                ('Gr0',  3.3e-4, True,   0.0,    1000,   None),
                ('E',    0,      True,   -1000,  1000,   None),
                ('v0',   43,     True,   -1e15,  1e15,   None),
                ('v1',   17.1,   True,   -1e15,  1e15,   None))
saveData(initParams, 'initParams')
from pyrho.datasets import *
ChR2data = loadChR2()
fitParams = fitModels(ChR2data, nStates=6, params=initParams, postFitOpt=True, relaxFact=2)
# Parameters automatically saved as 'fitted{}sParams'.format(nStates)

Example command-based NEURON simulation


In [ ]:
#setupNEURON(NEURONpath='/path/to/NEURON/') # Uncomment and set the path to compile NEURON
RhO = models['6']()
Prot = protocols['step']()
Prot.phis = [1e16, 1e15, 1e14]
Sim = simulators['NEURON'](Prot, RhO)
Sim.run()
Sim.plot()

Example custom stimulus protocol


In [ ]:
RhO = models['6']()
Prot = protocols['custom']()
Sim = simulators['Python'](Prot, RhO)

def pulseGenerator(run, phi, pulse):
    """Custom interpolation function for a step pulse with damped sinusoidal oscillations"""
    pStart, pEnd = pulse
    t = np.linspace(0, pEnd-pStart, 1001, endpoint=True)
    x = phi + 0.25*phi*np.sin(0.2*np.pi*t)*np.exp(-.05*t)
    return spline(pStart + t, x, ext=1, k=5)

Prot.phi_ft = pulseGenerator

Sim.run()
Sim.plot()

Example Brian simulation


In [ ]:
import brian2 as br

### Define the rhodopsin model
nStates = '6'
origParams = loadData('fitted{}sParams'.format(nStates))
#origParams = modelFits[nStates]['ChR2'] # Load from pre-fit models
RhO = models[nStates](origParams)

### Define the stimulation protocol
Prot = protocols['step']()
Prot.phis = [1e18]
Prot.Vs = [None]
Prot.cycles = [[150.,100.],[200.,75.]]

### Define network parameters
from brian2.units import Mohm
N = 80
pConnect = 0.2
psp = 'v_post += 1.5*mV'
delay = 3*ms
netParams = {'tau_m':10*ms, 'R_m':70*Mohm, 'E_m':-70*mV, 'v_t0':-50*mV, 'sigma':10*mV, 't_ref':4*ms}

### Define neuron model
eqRhO = '''dv/dt = ((-I*R_m)+E_m-v)/tau_m + sigma*xi*tau_m**-0.5 : volt''' + RhO.brian_phi_t
eqLIF = '''dv/dt = (E_m-v)/tau_m + sigma*xi*tau_m**-0.5 : volt'''

### Create neuron groups
G0 = br.NeuronGroup(N, eqRhO, threshold='v>v_t0', reset='v=E_m', refractory='t_ref', namespace=netParams, name='Inputs')
G0.v = 'rand()*(v_t0-E_m)+E_m'
G1 = br.NeuronGroup(N/2, eqLIF, threshold='v>v_t0', reset='v=E_m', refractory='t_ref', namespace=netParams)
G1.v = 'rand()*(v_t0-E_m)+E_m'
G2 = br.NeuronGroup(N/4, eqLIF, threshold='v>v_t0', reset='v=E_m', refractory='t_ref', namespace=netParams)
G2.v = 'rand()*(v_t0-E_m)+E_m'

### Create synapses
S1 = br.Synapses(G0, G1, pre=psp, delay=delay)
S1.connect(True, p=pConnect)
S2 = br.Synapses(G1, G2, pre=psp, delay=delay)
S2.connect(True, p=pConnect)

### Create monitors
monitors = {'states' : br.StateMonitor(G0, RhO.brianStateVars, record=0),   # Record states
            'I'      : br.StateMonitor(G0, 'I', record=0),                  # Record current
            'V'      : br.StateMonitor(G0, 'v', record=0),                  # Record voltage
            'spikes' : [br.SpikeMonitor(G0, name='Retina'),                 # Record spikes
                        br.SpikeMonitor(G1, name='LGN'), 
                        br.SpikeMonitor(G2, name='V1')]}

### Build the network
net = br.Network(br.collect())
net.add(monitors.values()) # #net.add(monitors)

### Run the simulation
Sim = simulators['Brian'](Prot, RhO, simParams['Brian'], net, netParams, monitors)
Sim.run()
Sim.plot()

In [ ]:
G0.equations  # Print out input layer equations

In [ ]: