diff --git a/pyrho/fitting.py b/pyrho/fitting.py index 38cd78c..d6ceb70 100755 --- a/pyrho/fitting.py +++ b/pyrho/fitting.py @@ -969,6 +969,249 @@ def solveGo(tlag, Gd, Go0=1000, tol=1e-9): +def fit6Kstates(fluxSet, quickSet, run, vInd, params, method=defMethod): # , verbose=config.verbose): + """ + fluxSet := ProtocolData set (of Photocurrent objects) to fit + quickSet:= ProtocolData set (of Photocurrent objects) with short pulses to fit opsin activation rates + run := Index for the run within the ProtocolData set + vInd := Index for Voltage clamp value within the ProtocolData set + params := Parameters object of model parameters with initial values [and bounds, expressions] + method := Fitting algorithm for the optimiser to use + """ + # verbose := Text output (verbosity) level + + + plotResult = bool(config.verbose > 1) + + nStates = '6K' + + ### Prepare the data + nRuns = fluxSet.nRuns + nPhis = fluxSet.nPhis + nVs = fluxSet.nVs + + assert(0 < nPhis) + assert(0 <= run < nRuns) + assert(0 <= vInd < nVs) + + Ions = [None for phiInd in range(nPhis)] + Ioffs = [None for phiInd in range(nPhis)] + tons = [None for phiInd in range(nPhis)] + toffs = [None for phiInd in range(nPhis)] + phis = [] + Is = [] + ts = [] + Vs = [] + + Icycles = [] + nfs = [] # Normalisation factors: e.g. /Ions[trial][-1] or /min(Ions[trial]) + + + # Trim off phase data + #frac = 1 + #chop = int(round(len(Ioffs[0])*frac)) + + for phiInd in range(nPhis): + targetPC = fluxSet.trials[run][phiInd][vInd] + #targetPC.alignToTime() + I = targetPC.I + t = targetPC.t + onInd = targetPC._idx_pulses_[0,0] ### Consider multiple pulse scenarios + offInd = targetPC._idx_pulses_[0,1] + Ions[phiInd] = I[onInd:offInd+1] + Ioffs[phiInd] = I[offInd:] #[I[offInd:] for I in Is] + tons[phiInd] = t[onInd:offInd+1]-t[onInd] + toffs[phiInd] = t[offInd:]-t[offInd] #[t[offInd:]-t[offInd] for t in ts] + #args=(Ioffs[phiInd][:chop+1],toffs[phiInd][:chop+1]) + phi = targetPC.phi + phis.append(phi) + + Is.append(I) + ts.append(t) + V = targetPC.V + Vs.append(V) + + Icycles.append(I[onInd:]) + nfs.append(I[offInd]) + #nfs.append(targetPC.I_peak_) + + + ### OFF PHASE + ### 3a. OFF CURVE: Fit biexponential to off curve to find lambdas + + OffKeys = ['Gd1', 'Gd2', 'Gf0', 'Ga3'] + + iOffPs = Parameters() # Create parameter dictionary + for k in OffKeys: + copyParam(k, params, iOffPs) + + ### Trim the first 10% of the off curve to allow I1 and I2 to empty? + + + ### This is an approximation based on the 4-state model which ignores the effects of Go1 and Go2 after light off. + + # lam1 + lam2 == Gd1 + Gd2 + Gf0 + Gb0 + # lam1 * lam2 == Gd1*Gd2 + Gd1*Gb0 + Gd2*Gf0 + + #, Gd2, Gf0, Gb0: (Gd1 + Gd2 + Gf0 + Gb0)/2 + #calcC = lambda b, Gd1, Gd2, Gf0, Gb0: np.sqrt(b**2 - (Gd1*Gd2 + Gd1*Gb0 + Gd2*Gf0)) + + def lams(p): + + Gd1 = p['Gd1'].value + Gd2 = p['Gd2'].value + Ga3 = p['Ga3'].value + + lam1 = Gd1 + lam2 = (Gd2 + Ga3) + return lam1, lam2 + + # Create dummy parameters for each phi + for phiInd in range(nPhis): + Iss = Ioffs[phiInd][0] + if Iss < 0: + iOffPs.add('Islow_'+str(phiInd), value=0.2*Iss, vary=True, max=0) + iOffPs.add('Ifast_'+str(phiInd), value=0.8*Iss, vary=True, max=0, expr='{} - {}'.format(Iss, 'Islow_'+str(phiInd))) + else: + iOffPs.add('Islow_'+str(phiInd), value=0.2*Iss, vary=True, min=0) + iOffPs.add('Ifast_'+str(phiInd), value=0.8*Iss, vary=True, min=0, expr='{} - {}'.format(Iss, 'Islow_'+str(phiInd))) + + def fit6Koff(p,t,trial): + Islow = p['Islow_'+str(trial)].value + Ifast = p['Ifast_'+str(trial)].value + lam1, lam2 = lams(p) + return Islow*np.exp(-lam1*t) + Ifast*np.exp(-lam2*t) + + def err6Koff(p,Ioffs,toffs): + """Normalise by the first element of the off-curve""" # [-1] + return np.r_[ [(Ioffs[i] - fit6Koff(p,toffs[i],i))/Ioffs[i][0] for i in range(len(Ioffs))] ] + + #fitfunc = lambda p, t: -(p['a0'].value + p['a1'].value*np.exp(-lams(p)[0]*t) + p['a2'].value*np.exp(-lams(p)[1]*t)) + ##fitfunc = lambda p, t: -(p['a0'].value + p['a1'].value*np.exp(-p['lam1'].value*t) + p['a2'].value*np.exp(-p['lam2'].value*t)) + #errfunc = lambda p, Ioff, toff: Ioff - fitfunc(p,toff) + + offPmin = minimize(err6Koff, iOffPs, args=(Ioffs,toffs), method=method)#, fit_kws={'maxfun':100000}) + pOffs = offPmin.params + + reportFit(offPmin, "Off-phase fit report for the 6K-state model", method) + if config.verbose > 0: + print('Gd1 = {}; Gd2 = {}; Gf0 = {}'.format(pOffs['Gd1'].value, pOffs['Gd2'].value, + pOffs['Gf0'].value)) + + if plotResult: + lam1, lam2 = lams(pOffs) + plotOffPhaseFits(toffs, Ioffs, pOffs, phis, nStates, fit6Koff, lam1, lam2, Gd=None) + + + # Fix off-curve parameters + for k in OffKeys: + pOffs[k].vary = False + + + ### Calculate Go (1/tau_opsin) + print('\nCalculating opsin activation rate') + # Assume that Gd1 > Gd2 + # Assume that Gd = Gd1 for short pulses + + def solveGo(tlag, Gd, Go0=1000, tol=1e-9): + Go, Go_m1 = Go0, 0 + while abs(Go_m1 - Go) > tol: + Go_m1 = Go + Go = ((tlag*Gd) - np.log(Gd/Go_m1))/tlag + #Go_m1, Go = Go, ((tlag*Gd) - np.log(Gd/Go_m1))/tlag + return Go + + #if 'shortPulse' in dataSet: # Fit Go + if quickSet.nRuns > 1: + #from scipy.optimize import curve_fit + # Fit tpeak = tpulse + tmaxatp0 * np.exp(-k*tpulse) + #dataSet['shortPulse'].getProtPeaks() + #tpeaks = dataSet['shortPulse'].IrunPeaks + + #PD = dataSet['shortPulse'] + PCs = [quickSet.trials[p][0][0] for p in range(quickSet.nRuns)] # Aligned to the pulse i.e. t_on = 0 + #[pc.alignToTime() for pc in PCs] + + #tpeaks = np.asarray([PD.trials[p][0][0].tpeak for p in range(PD.nRuns)]) # - PD.trials[p][0][0].t[0] + #tpulses = np.asarray([PD.trials[p][0][0].Dt_ons[0] for p in range(PD.nRuns)]) + tpeaks = np.asarray([pc.t_peak_ for pc in PCs]) + tpulses = np.asarray([pc.Dt_ons_[0] for pc in PCs]) + + devFunc = lambda tpulses, t0, k: tpulses + t0 * np.exp(-k*tpulses) + p0 = (0, 1) + popt, pcov = curve_fit(devFunc, tpulses, tpeaks, p0=p0) + if plotResult: + fig = plt.figure() + ax = fig.add_subplot(111, aspect='equal') + nPoints = 10*int(round(max(tpulses))+1) # 101 + tsmooth = np.linspace(0, max(tpulses), nPoints) + ax.plot(tpulses, tpeaks, 'x') + ax.plot(tsmooth, devFunc(tsmooth, *popt)) + ax.plot(tsmooth, tsmooth, '--') + ax.set_ylim([0, max(tpulses)]) #+5 + ax.set_xlim([0, max(tpulses)]) #+5 + #plt.tight_layout() + #plt.axis('equal') + plt.show() + + # Solve iteratively Go = ((tlag*Gd) - np.log(Gd/Go))/tlag + Gd1 = pOffs['Gd1'].value + Go = solveGo(tlag=popt[0], Gd=Gd1, Go0=1000, tol=1e-9) + print('t_lag = {:.3g}; Gd = {:.3g} --> Go = {:.3g}'.format(popt[0], Gd1, Go)) + + elif quickSet.nRuns == 1: #'delta' in dataSet: + #PD = dataSet['delta'] + #PCs = [PD.trials[p][0][0] for p in range(PD.nRuns)] + PC = quickSet.trials[0][0][0] + tlag = PC.Dt_lag_ # := Dt_lags_[0] ############################### Add to Photocurrent... + Go = solveGo(tlag=tlag, Gd=Gd1, Go0=1000, tol=1e-9) + print('t_lag = {:.3g}; Gd = {:.3g} --> Go = {:.3g}'.format(tlag, Gd1, Go)) + + else: + Go = 1 # Default + print('No data found to estimate Go: defaulting to Go = {}'.format(Go)) + + + ### ON PHASE + + iOnPs = Parameters() # deepcopy(params) + + # Set parameters from Off-curve optimisation + for k in OffKeys: + copyParam(k, pOffs, iOnPs) + + # Set parameters from general rhodopsin analysis routines + for k in ['Go1', 'Go2', 'k1', 'k2', 'k3', 'k_f', 'k_b', 'gam', 'p', 'q', 'phi_m', 'g0', 'Gb', 'E', 'v0', 'v1']: #.extend(OffKeys): + copyParam(k, params, iOnPs) + + # Set parameters from short pulse calculations + iOnPs['Go1'].value = Go; iOnPs['Go1'].vary = False + iOnPs['Go2'].value = Go; iOnPs['Go2'].vary = False + + RhO = models['6K']() + + ### Trim down ton? Take 10% of data or one point every ms? ==> [0::5] + + if config.verbose > 2: + print('Optimising ',end='') + + onPmin = minimize(errOnPhase, iOnPs, args=(Ions,tons,RhO,Vs,phis), method=method) + pOns = onPmin.params + + reportFit(onPmin, "On-phase fit report for the 6K-state model", method) + + if config.verbose > 0: + print('k1 = {}; k2 = {}; k_f = {}; k_b = {}'.format(pOns['k1'].value, pOns['k2'].value, + pOns['k_f'].value, pOns['k_b'].value)) + print('gam = {}; phi_m = {}; p = {}; q = {}'.format(pOns['gam'].value, pOns['phi_m'].value, + pOns['p'].value, pOns['q'].value)) + + fitParams = pOns + + return fitParams, onPmin + + + #TODO: Tidy up and refactor getRecoveryPeaks and fitRecovery def getRecoveryPeaks(recData, phiInd=None, vInd=None, usePeakTime=False): @@ -1769,11 +2012,14 @@ def fitModel(dataSet, nStates='3', params=None, postFitOpt=True, relaxFact=2, me """Fit a model (with initial parameters) to a dataset of optogenetic photocurrents.""" ### Define non-optimised parameters to exclude in post-fit optimisation - nonOptParams = ['Gr0', 'E', 'v0', 'v1'] - if not isinstance(nStates, str): nStates = str(nStates) # .lower() + if nStates == '3' or nStates == '4' or nStates == '6': + nonOptParams = ['Gr0', 'E', 'v0', 'v1'] + elif nStates == '6K': + nonOptParams = ['E', 'v0', 'v1', 'Ga3'] + if nStates not in modelParams: print(f"Error in selecting model {nStates} - please choose from {list(modelParams)} states") raise NotImplementedError(nStates) @@ -2000,6 +2246,11 @@ def fitModel(dataSet, nStates='3', params=None, postFitOpt=True, relaxFact=2, me constrainedParams = ['Gd1', 'Gd2', 'Gf0', 'Gb0', 'Go1', 'Go2'] #constrainedParams = ['Go1', 'Go2', 'Gf0', 'Gb0'] #nonOptParams.append(['Gd1', 'Gd2']) + elif nStates == '6K': + fittedParams, miniObj = fit6Kstates(setPC, quickSet, runInd, vIndm70, fitParams, method) # , verbose) + constrainedParams = ['Gd1', 'Gd2', 'Gf0', 'Go1', 'Go2', 'Gb'] + #constrainedParams = ['Go1', 'Go2', 'Gf0', 'Gb0'] + #nonOptParams.append(['Gd1', 'Gd2']) else: raise Exception(f'Invalid choice for nStates: {nStates}!') diff --git a/pyrho/models.py b/pyrho/models.py index 52675dd..4be08f5 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -36,7 +36,10 @@ class RhodopsinModel(PyRhOobject): def __init__(self, params=None, rhoType=rhoType): if params is None: - params = modelParams[str(self.nStates)] + try: + params = modelParams[str(self.nStates)] + except: + pass self.rhoType = rhoType # E.g. 'ChR2' or 'ArchT' self.setParams(params) @@ -89,6 +92,12 @@ def initStates(self, phi, s0=None): """Clear state arrays and set transition rates.""" if s0 is None: s0 = self.s_0 + #assert(len(s0) == self.nStates) + #if self.nStates=='6K': + # self.states = np.vstack((np.empty([0, 6]), s0)) + # + #else: + # self.states = np.vstack((np.empty([0, self.nStates]), s0)) assert len(s0) == self.nStates self.states = np.vstack((np.empty([0, self.nStates]), s0)) self.t = [0] @@ -819,11 +828,191 @@ def calcfphi(self, states=None): def calcSoln(self, t, s0=None): raise NotImplementedError(self.nStates) +class RhO_6Kstates(RhodopsinModel): + """Class definition for the 6K-state model""" + + # Class attributes + nStates = 6 + useAnalyticSoln = False + s_0 = np.array([1, 0, 0, 0, 0, 0]) # [s1_0=1, s2_0=0, s3_0=0, s4_0=0, s5_0=0, s6_0=0] # array not necessary + phi_0 = 0.0 # Default initial flux + stateVars = ['C1', 'I1', 'O1', 'O2', 'I2', 'C2'] # stateVars[0] is the 'ground' state + stateLabels = ['$C_1$', '$I_1$', '$O_1$', '$O_2$', '$I_2$', '$C_2$'] + photoFuncs = ['_calcGa1', '_calcGa2', '_calcGf'] + photoRates = ['Ga1', 'Ga2', 'Gf'] + photoLabels = ['$G_{a1}$', '$G_{a2}$', '$G_{f}$'] + constRates = ['Go1', 'Go2', 'Ga3', 'Gd1', 'Gd2', 'Gb'] + constLabels = ['$G_{o1}$', '$G_{o2}$', '$G_{a3}$', '$G_{d1}$', '$G_{d2}$', '$G_{b}$'] + + paramsList = ['g0', 'gam', 'phi_m', 'k1', 'k2', 'k3', 'p', + 'Gf0', 'k_f', 'Gb0', 'k_b', 'q', 'Go1', 'Go2', + 'Ga3', 'Gd1', 'Gd2', 'Gb', 'E', 'v0', 'v1'] # List of model constants + + connect = [[0, 0, 1, 1, 0, 1], # s_1 --> s_i=1...6 + [1, 0, 1, 0, 0, 0], # s_2 --> + [0, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1], + [1, 0, 0, 1, 0, 0]] + + equations = r""" + $$ \dot{C_1} = G_{d1}O_1 + G_{b}C_2 + G_{a3}O_2 - (G_{a1} + G_{f})(\phi)C_1 $$ + $$ \dot{I_1} = G_{a1}(\phi)C_1 - G_{o1}I_1 $$ + $$ \dot{O_1} = G_{o1}I_1 - G_{d1}O_1 $$ + $$ \dot{O_2} = G_{o2}I_2 - (G_{d2} + G_{a3})O_2 $$ + $$ \dot{I_2} = G_{a2}(\phi)C_2 - G_{o2}I_2 $$ + $$ \dot{C_2} = G_{d2}O_2 + G_{f}(\phi)C1 - (G_{b} + G_{a2}(\phi))C_2 $$ + $$ C_1 + I_1 + O_1 + O_2 + I_2 + C_2 = 1 $$ + $$$$ + $$ G_{a1}(\phi) = k_{1} \frac{\phi^p}{\phi^p + \phi_m^p} $$ + $$ G_{a2}(\phi) = k_{2} \frac{\phi^p}{\phi^p + \phi_m^p} $$ + $$ G_{a3}(\phi) = k_{3} \frac{\phi^p}{\phi^p + \phi_m^p} $$ + $$ G_{f}(\phi) = k_{f} \frac{\phi^q}{\phi^q + \phi_m^q} + G_{f0} $$ + $$ G_{b}(\phi) = k_{b} \frac{\phi^q}{\phi^q + \phi_m^q} + G_{b0} $$ + $$$$ + $$ f_{\phi}(\phi) = O_1+\gamma O_2 $$ + $$ f_v(v) = v_1\frac{1-e^{-(v-E)/v_0}}{(v-E)} $$ + $$ I_{\phi} = g_0 \cdot f_{\phi}(\phi) \cdot f_v(v) \cdot (v-E) $$ + """ + + brianStateVars = ['C_1', 'I_1', 'O_1', 'O_2', 'I_2', 'C_2'] + + brian = ''' + dC_1/dt = Gd1*O_1 + Gb*C_2 + Ga3*O_2 - (Ga1+Gf)*C_1 : 1 + dI_1/dt = Ga1*C_1 - Go1*I_1 : 1 + dO_1/dt = Go1*I_1 - Gd1*O_1 : 1 + dO_2/dt = Go2*I_2 - (Gd2+Ga3)*O_2 : 1 + dI_2/dt = Ga2*C_2 - Go2*I_2 : 1 + C_2 = 1 - C_1 - I_1 - O_1 - O_2 - I_2 : 1 + H_p = Theta*((phi**p)/(phi**p+phi_m**p)) : 1 + H_q = Theta*((phi**q)/(phi**q+phi_m**q)) : 1 + Ga1 = k1*H_p : second**-1 + Ga2 = k2*H_p : second**-1 + Ga3 = k3*H_p : second**-1 + Gf = k_f*H_q + Gf0 : second**-1 + Gb = k_b*H_q + Gb0 : second**-1 + f_v = (1-exp(-(v-E)/v0))/((v-E)/v1) : 1 + f_phi = O_1+gam*O_2 : 1 + I = g0*f_phi*f_v*(v-E) : amp + phi : metre**-2*second**-1 (shared) + Theta = int(phi > 0*phi) : 1 (shared) + ''' + + brian_phi_t = ''' + dC_1/dt = Gd1*O_1 + Gb*C_2 + Ga3*O_2 - (Ga1+Gf)*C_1 : 1 + dI_1/dt = Ga1*C_1 - Go1*I_1 : 1 + dO_1/dt = Go1*I_1 - Gd1*O_1 : 1 + dO_2/dt = Go2*I_2 - (Gd2+Ga3)*O_2 : 1 + dI_2/dt = Ga2*C_2 - Go2*I_2 : 1 + C_2 = 1 - C_1 - I_1 - O_1 - O_2 - I_2 : 1 + Theta = int(phi(t) > 0*phi(t)) : 1 (shared) + H_p = Theta*((phi(t)**p)/(phi(t)**p+phi_m**p)) : 1 + H_q = Theta*((phi(t)**q)/(phi(t)**q+phi_m**q)) : 1 + Ga1 = k1*H_p : second**-1 + Ga2 = k2*H_p : second**-1 + Ga3 = k3*H_p : second**-1 + Gf = k_f*H_q + Gf0 : second**-1 + Gb = k_b*H_q + Gb0 : second**-1 + f_v = (1-exp(-(v-E)/v0))/((v-E)/v1) : 1 + f_phi = O_1+gam*O_2 : 1 + I = g0*f_phi*f_v*(v-E) : amp + ''' + + def _calcGa1(self, phi): + #return self.a10*(phi/self.phi0) + return self.k1 * phi**self.p/(phi**self.p + self.phi_m**self.p) + + def _calcGa2(self, phi): + #return self.b40*(phi/self.phi0) + return self.k2 * phi**self.p/(phi**self.p + self.phi_m**self.p) + + def _calcGf(self, phi): + #return self.a30 + self.a31*np.log(1+(phi/self.phi0)) + return self.Gf0 + self.k_f * phi**self.q/(phi**self.q + self.phi_m**self.q) + + def setLight(self, phi): + if phi < 0: + phi = 0 + self.phi = phi + self.Ga1 = self._calcGa1(phi) + self.Ga2 = self._calcGa2(phi) + self.Gf = self._calcGf(phi) + if config.verbose > 1: + self.dispRates() + + def dispRates(self): # This needs rechecking + """Print photosensitive transition rates""" + print("Transition rates (phi={:.3g}): C1 --[Ga1={:.3g}]--> O1 --[Gf={:.3g}]--> O2".format(self.phi, self.Ga1, self.Gf)) + #print("Transition rates (phi={:.3g}): O1 <--[Gb={:.3g}]-- O2 <--[Ga2={:.3g}]-- C2".format(self.phi, self.Gb, self.Ga2)) + + def solveStates(self, s_0, t, phi_t=None): + """Differential equations of the 6-state model to be solved by odeint""" + if phi_t is not None: + self.setLight(float(phi_t(t))) + C1, I1, O1, O2, I2, C2 = s_0 # Unpack state vector + dC1dt = -(self.Ga1+self.Gf)*C1 + self.Gd1*O1 + self.Ga3*O2 + self.Gb*C2 + dI1dt = self.Ga1*C1 - self.Go1*I1 + dO1dt = self.Go1*I1 - self.Gd1*O1 + dO2dt = -(self.Ga3+self.Gd2)*O2 + self.Go2*I2 + dI2dt = -self.Go2*I2 + self.Ga2*C2 + dC2dt = self.Gf*C1 + self.Gd2*O2 - (self.Ga2+self.Gb)*C2 + return np.array([dC1dt, dI1dt, dO1dt, dO2dt, dI2dt, dC2dt]) + + def jacobian(self, s_0, t, phi_t=None): + return np.array([[-(self.Ga1+self.Gf), 0, self.Gd1, self.Ga3, 0, self.Gb], + [self.Ga1, -self.Go1, 0, 0, 0, 0], + [0, self.Go1, -self.Gd1, 0, 0, 0], + [0, 0, 0, -(self.Ga3+self.Gd2), self.Go2, 0], + [0, 0, 0, 0, -self.Go2, self.Ga2], + [self.Gf, 0, 0, self.Gd2, 0, -(self.Ga2+self.Gb)]]) + + def hessian(self, s_0, t): + """ + Hessian matrix for scipy.optimize.minimize. + (Only for Newton-CG, dogleg, trust-ncg.) + H(f)_ij(X) = D_iD_jf(X) + """ + return np.zeros((6, 6)) + + def calcSteadyState(self, phi): + self.setLight(phi) + Ga1 = self.Ga1 + Go1 = self.Go1 + Gf = self.Gf + Gd2 = self.Gd2 + Gd1 = self.Gd1 + Gb = self.Gb + Go2 = self.Go2 + Ga2 = self.Ga2 + Ga3 = self.Ga3 + + denom = (Ga1*Ga2*Ga3*Gd1*Go2 + Ga1*Ga2*Ga3*Go1*Go2 + Ga1*Ga3*Gb*Gd1*Go2 + Ga1*Ga3*Gb*Go1*Go2 + Ga1*Gb*Gd1*Gd2*Go2 + Ga1*Gb*Gd2*Go1*Go2 + Ga2*Ga3*Gd1*Gf*Go1 + Ga2*Ga3*Gd1*Go1*Go2 + Ga2*Gd1*Gd2*Gf*Go1 + Ga2*Gd1*Gf*Go1*Go2 + Ga3*Gb*Gd1*Go1*Go2 + Ga3*Gd1*Gf*Go1*Go2 + Gb*Gd1*Gd2*Go1*Go2 + Gd1*Gd2*Gf*Go1*Go2) + C1ss = (Ga2*Ga3*Gd1*Go1*Go2 + Ga3*Gb*Gd1*Go1*Go2 + Gb*Gd1*Gd2*Go1*Go2) + I1ss = (Ga1*Ga2*Ga3*Gd1*Go2 + Ga1*Ga3*Gb*Gd1*Go2 + Ga1*Gb*Gd1*Gd2*Go2) + O1ss = (Ga1*Ga2*Ga3*Go1*Go2 + Ga1*Ga3*Gb*Go1*Go2 + Ga1*Gb*Gd2*Go1*Go2) + O2ss = (Ga2*Gd1*Gf*Go1*Go2) + I2ss = (Ga2*Ga3*Gd1*Gf*Go1 + Ga2*Gd1*Gd2*Gf*Go1) + C2ss = (Ga3*Gd1*Gf*Go1*Go2 + Gd1*Gd2*Gf*Go1*Go2) + self.steadyStates = np.array([C1ss, I1ss, O1ss, O2ss, I2ss, C2ss]) / denom + return self.steadyStates + + def calcfphi(self, states=None): + """Calculate the conductance scalar from the photocycle""" + if states is None: + states = self.states + C1, I1, O1, O2, I2, C2 = states.T + gam = self.gam + return O1 + gam * O2 + + def calcSoln(self, t, s0=None): + raise NotImplementedError(self.nStates) + models = { '3': RhO_3states, 3: RhO_3states, '4': RhO_4states, 4: RhO_4states, - '6': RhO_6states, 6: RhO_6states + '6': RhO_6states, 6: RhO_6states, + '6K': RhO_6Kstates } @@ -835,6 +1024,8 @@ def selectModel(nStates): return RhO_4states() elif int(nStates) == 6 or nStates == 'six': return RhO_6states() + elif nStates == '6K': + return RhO_6Kstates() else: print("Error in selecting model - please choose from 3, 4 or 6 states") raise NotImplementedError(nStates) diff --git a/pyrho/parameters.py b/pyrho/parameters.py index 0fdaa91..56cc679 100755 --- a/pyrho/parameters.py +++ b/pyrho/parameters.py @@ -75,18 +75,20 @@ # Default model parameters -modelParams = {'3': Parameters(), '4': Parameters(), '6': Parameters()} +modelParams = {'3': Parameters(), '4': Parameters(), '6': Parameters(), '6K': Parameters()} modelList = list(modelParams) # List of keys: list(modelParams.keys()) #This could be removed stateLabs = {'3': 'Three', 3: 'Three', '4': 'Four', 4: 'Four', - '6': 'Six', 6: 'Six'} + '6': 'Six', 6: 'Six', + '6K': 'SixK'} modelFits = { '3': {'ChR2': Parameters(), 'NpHR': Parameters(), 'ArchT': Parameters()}, '4': {'ChR2': Parameters()}, - '6': {'ChR2': Parameters()} + '6': {'ChR2': Parameters()}, + '6K': {'ChR2': Parameters()} } # Replace with defaultdict with default=key @@ -217,12 +219,35 @@ ('v0', 43, True, -1e15, 1e15,None), ('v1', 17.1, True, -1e15, 1e15,None)) +modelFits['6K']['ChR2'].add_many( + ('g0', 2.52e4, True, 0.0, 1e15, None), + ('gam', 0.0161, True, 0.0, 1, None), + ('phi_m', 3.54e17,True, 1e15, 1e19, None), + ('k1', 13.4, True, 0.0, 1000, None), + ('k2', 2.71, True, 0.0, 1000, None), + ('k3', 2.71, True, 0.0, 1000, None), + ('p', 0.985, True, 0.1, 5, None), + ('Gf0', 0.0389, True, 0.0, 1000, None), + ('k_f', 0.103, True, 0.0, 1000, None), + ('k_b', 0.139, True, 0.0, 1000, None), + ('q', 1.58, True, 0.1, 5, None), + ('Go1', 2, True, 0.0, 1000, None), + ('Go2', 0.0567, True, 0.0, 1000, None), + ('Gd1', 0.112, True, 0.0, 1000, None), + ('Gd2', 0.0185, True, 0.0, 1000, None), + ('Ga3', 250, True, 0.0, 500, None), + ('Gb', 40000, True, 3400.0, 45000, None), + ('E', 0, True, -1000,1000, None), + ('v0', 43, True, -1e15, 1e15,None), + ('v1', 17.1, True, -1e15, 1e15,None)) + defaultOpsinType = 'ChR2' rhoType = defaultOpsinType # Set this when selecting modelParams['3'] = modelFits['3'][defaultOpsinType] modelParams['4'] = modelFits['4'][defaultOpsinType] modelParams['6'] = modelFits['6'][defaultOpsinType] +modelParams['6K'] = modelFits['6K'][defaultOpsinType] unitPrefixes = {} # Use a units library to convert between different prefixes diff --git a/pyrho/run_sims.py b/pyrho/run_sims.py new file mode 100644 index 0000000..eb831bd --- /dev/null +++ b/pyrho/run_sims.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jun 21 23:34:21 2022 + +@author: beers +""" + +from pyrho import * +from pyrho.datasets import * + +initParams = Parameters() +initParams.add_many( + # Name Value Vary Min Max Expr + ('g0', 2.52e4, True, 0.0, 1e15, None), + ('gam', 0.0161, True, 0.0, 1, None), # Max=1 if gO1 >= gO2 + ('phi_m', 3.54e17,True, 1e15, 1e19, None), + ('k1', 13.4, True, 0.0, 1000, None), + ('k2', 2.71, True, 0.0, 1000, None), + ('k3', 2.71, True, 0.0, 1000, None), # New: find values + ('p', 0.985, True, 0.1, 5, None), + ('Gf0', 0.0389, True, 0.0, 1000, None), + ('k_f', 0.103, True, 0.0, 1000, None), + #('Gb0', 0.0198, True, 0.0, 1000, None), + ('k_b', 0.139, True, 0.0, 1000, None), + ('q', 1.58, True, 0.1, 5, None), + ('Go1', 2, True, 0.0, 1000, None), + ('Go2', 0.0567, True, 0.0, 1000, None), + ('Gd1', 0.112, True, 0.0, 1000, None), + ('Gd2', 0.0185, True, 0.0, 1000, None), + ('Ga3', 250, True, 0.0, 500, None), + ('Gb', 40000, True, 3400.0, 45000, None), + ('E', 0, True, -1000,1000, None), + ('v0', 43, True, -1e15, 1e15,None), + ('v1', 17.1, True, -1e15, 1e15,None)) + +saveData(initParams, 'initParams') + +ChR2data = loadChR2() +fitParams = fitModels(ChR2data, nStates='6K', params=initParams, postFitOpt=True, relaxFact=2) \ No newline at end of file