From a1ebbb2539902ba1a6e4970ea74de34eaded73b1 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:01:59 +0100 Subject: [PATCH 01/12] Include naming exception in model base class Include exception in RhodopsinModel base class when nStates='6K' and cannot be converted to string --- pyrho/models.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyrho/models.py b/pyrho/models.py index df8a59b..b3d3da6 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -38,7 +38,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) From c95278b42bc882affe651ff76cc83e3f56db91db Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:11:59 +0100 Subject: [PATCH 02/12] Update models.py Remove assertion that 'nStates' is an int, and account for this when initialising 'states' --- pyrho/models.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pyrho/models.py b/pyrho/models.py index b3d3da6..51f130d 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -93,8 +93,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) - self.states = np.vstack((np.empty([0, self.nStates]), s0)) + #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)) self.t = [0] self.pulseInd = np.empty([0, 2], dtype=int) # Light on and off indexes for each pulse self.ssInf = [] From 29ac814d24e3e71698b88a864ae19f7147b5509d Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:14:28 +0100 Subject: [PATCH 03/12] Create a Kuhne model class --- pyrho/models.py | 179 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 179 insertions(+) diff --git a/pyrho/models.py b/pyrho/models.py index 51f130d..ec15403 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -990,6 +990,185 @@ 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 = '6K' + 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 = OrderedDict([('3', RhO_3states), (3, RhO_3states), ('4', RhO_4states), (4, RhO_4states), From 6488cdc16ec450592f2439e74191331b592c0e72 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:17:41 +0100 Subject: [PATCH 04/12] Update models.py Include the 6K model in the ordered dict containing all models --- pyrho/models.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyrho/models.py b/pyrho/models.py index ec15403..d3049d5 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -1172,7 +1172,8 @@ def calcSoln(self, t, s0=None): models = OrderedDict([('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)]) def selectModel(nStates): @@ -1183,6 +1184,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) From 353d6d76bb6afc6cbebfd0dce0fbf7d386130525 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:24:08 +0100 Subject: [PATCH 05/12] Update parameters.py Include Kuhne model within the list of default model parameters --- pyrho/parameters.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pyrho/parameters.py b/pyrho/parameters.py index 1de263e..350d013 100755 --- a/pyrho/parameters.py +++ b/pyrho/parameters.py @@ -77,17 +77,19 @@ # Default model parameters -modelParams = OrderedDict([('3', Parameters()), ('4', Parameters()), ('6', Parameters())]) +modelParams = OrderedDict([('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 = OrderedDict([('3', OrderedDict([('ChR2', Parameters()), ('NpHR', Parameters()), ('ArchT', Parameters())])), ('4', OrderedDict([('ChR2', Parameters())])), - ('6', OrderedDict([('ChR2', Parameters())]))]) + ('6', OrderedDict([('ChR2', Parameters())])), + ('6K', OrderedDict([('ChR2', Parameters())]))]) # Replace with defaultdict with default=key modelLabels = OrderedDict([('E', 'E'), ('g0', 'g_0'), ('p', 'p'), From c5984f54dd7822e8496a5b3a42cc56118041295b Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:25:44 +0100 Subject: [PATCH 06/12] Update parameters.py Define default parameters for 6K model --- pyrho/parameters.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/pyrho/parameters.py b/pyrho/parameters.py index 350d013..218d459 100755 --- a/pyrho/parameters.py +++ b/pyrho/parameters.py @@ -214,6 +214,28 @@ ('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 From 22003e688e79d6fdb4ea922ca8c212d9b307badc Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:27:38 +0100 Subject: [PATCH 07/12] Update parameters.py Include 6K model in 'modelParams' --- pyrho/parameters.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyrho/parameters.py b/pyrho/parameters.py index 218d459..ca93093 100755 --- a/pyrho/parameters.py +++ b/pyrho/parameters.py @@ -242,6 +242,7 @@ 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 From 1dbb410bf6cf6965970112467ec582dda846567c Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:30:57 +0100 Subject: [PATCH 08/12] Update fitting.py Describe fitting procedure for 6K model --- pyrho/fitting.py | 242 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 242 insertions(+) diff --git a/pyrho/fitting.py b/pyrho/fitting.py index 5bd161d..e170f2e 100755 --- a/pyrho/fitting.py +++ b/pyrho/fitting.py @@ -963,6 +963,248 @@ 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): From 24fa233c75a893f0943f50b66b574e7437a0d2f1 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:38:44 +0100 Subject: [PATCH 09/12] Update fitting.py Redefine list of non-optimised parameters for 6K model. Warning: This removes the check to ensure only one of the possible models (3, 4, 6) has been selected. To do: Ensure only one of the four models is selected. --- pyrho/fitting.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/pyrho/fitting.py b/pyrho/fitting.py index e170f2e..4e17b13 100755 --- a/pyrho/fitting.py +++ b/pyrho/fitting.py @@ -1206,6 +1206,7 @@ def solveGo(tlag, Gd, Go0=1000, tol=1e-9): + #TODO: Tidy up and refactor getRecoveryPeaks and fitRecovery def getRecoveryPeaks(recData, phiInd=None, vInd=None, usePeakTime=False): """Aggregate the times and currents of the photocurrent peaks for fitting Gr0 @@ -2007,12 +2008,14 @@ def fitModel(dataSet, nStates=3, params=None, postFitOpt=True, relaxFact=2, meth ### Define non-optimised parameters to exclude in post-fit optimisation - nonOptParams = ['Gr0', 'E', 'v0', 'v1'] - - if isinstance(nStates, str): - nStates = int(nStates) # .lower() - if nStates == 3 or nStates == 4 or nStates == 6: + try: + nStates = int(nStates) + except: pass + if nStates == 3 or nStates == 4 or nStates == 6: + nonOptParams = ['Gr0', 'E', 'v0', 'v1'] + elif nStates == '6K': + nonOptParams = ['E', 'v0', 'v1', 'Ga3'] else: print("Error in selecting model - please choose from 3, 4 or 6 states") raise NotImplementedError(nStates) From 545386541216ea01de66c97e4961704c77fb30d6 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 00:49:37 +0100 Subject: [PATCH 10/12] Create run_sims.py --- pyrho/run_sims.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 pyrho/run_sims.py 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 From 9e1366c3b6d9d7b5f523088c71dddef724bb23f1 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Wed, 22 Jun 2022 16:21:41 +0100 Subject: [PATCH 11/12] Update fitting.py Include constrained parameters for 6K model. --- pyrho/fitting.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pyrho/fitting.py b/pyrho/fitting.py index 4e17b13..f8c660e 100755 --- a/pyrho/fitting.py +++ b/pyrho/fitting.py @@ -2242,6 +2242,11 @@ def fitModel(dataSet, nStates=3, params=None, postFitOpt=True, relaxFact=2, meth 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('Invalid choice for nStates: {}!'.format(nStates)) From e99140a22cdf50087a72e74cb84356fd75de2346 Mon Sep 17 00:00:00 2001 From: AdrianBeersingVasquez <57265495+AdrianBeersingVasquez@users.noreply.github.com> Date: Tue, 18 Apr 2023 11:16:20 +0100 Subject: [PATCH 12/12] Tidy up comment --- pyrho/models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyrho/models.py b/pyrho/models.py index 367fd51..4be08f5 100755 --- a/pyrho/models.py +++ b/pyrho/models.py @@ -832,7 +832,7 @@ class RhO_6Kstates(RhodopsinModel): """Class definition for the 6K-state model""" # Class attributes - nStates = 6 # '6K' + 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