Code for generating 24_10_24 protocol

Author: Frankie Patten-Elliott

methods/classes.py

@@ -59,7 +59,7 @@ def __init__(self, model, protocol, times, win, conc, param_dict):
         self._conc = conc
         # TODO currently hardcoded to get number of pulses
         if model.split("-")[0] == 'sis':
-            self.n_pulses = 10
+            self.n_pulses = 20
         elif times[-1] != 14999.5:
             self.n_pulses = int(np.floor(250000/times[-1]))
         else:

methods/funcs.py

@@ -18,8 +18,10 @@ def generate_data(herg_model, drug_vals, prot, sd, max_time, bounds, m_sel, conc
     # TODO currently hardcoded to get number of sweeps
     if max_time != 15e3 and herg_model != '2024_Joey_sis_25C':
         swps = int(np.floor(250000/max_time))
-    else:
+    elif max_time != 15350:
         swps = sweeps
+    else:
+        swps = 20
 
     # define protocol
     protocol = myokit.load_protocol(prot)
@@ -161,11 +163,12 @@ def model_outputs(model_pars, herg_model, prot, times, concs, alt_protocol = Non
     # get herg parameters
     if herg_model == '2019_37C':
         herg_vals = [2.07e-3, 7.17e-2, 3.44e-5, 6.18e-2, 4.18e-1, 2.58e-2, 4.75e-2, 2.51e-2, 33.3]
-    elif herg_model == 'kemp':
+    elif herg_model == 'kemp' or herg_model == '2024_Joey_sis_25C':
         herg_vals = []
     model_out = {}
     # get no. of sweeps such that the total length is approximately 250s
-    swps = int(np.floor(250000/(times[-1]+alt_times[-1])))
+    #swps = int(np.floor(250000/(times[-1]+alt_times[-1])))
+    swps = 5
     win = (times >= wins[0]) & (times < wins[1])
     if alt_protocol is not None:
         win_alt = (alt_times >= wins_alt[0]) & (alt_times < wins_alt[1])
@@ -181,39 +184,44 @@ def model_outputs(model_pars, herg_model, prot, times, concs, alt_protocol = Non
                                 analytical=True)
                 # set hERG model parameters
                 model.set_fix_parameters({p:v for p, v in zip(herg_pars, herg_vals)})
-            else:
+            elif herg_model == 'kemp':
                 model = Model(f'kemp-m{m}',
                                 prot,
                                 parameters=['binding'],
                                 analytical=True)
+            elif herg_model == '2024_Joey_sis_25C':
+                model = Model(f'sis-m{m}',
+                                prot,
+                                parameters=['binding'],
+                                analytical=True)
             # fix kt if necessary
             if m in ['12', '13']:
                 model.fix_kt()
             try:
                 # loop to simulate and append proportion open model output for no. of sweeps 
                 model_milnes = []
-                control = []
-                model.set_dose(0)
-                before = model.simulate(binding_params, times)[win]
-                before_alt = model.simulate(binding_params, alt_times)[win_alt]
+                #control = []
+                #model.set_dose(0)
+                #before = model.simulate(binding_params, times)[win]
+                #before_alt = model.simulate(binding_params, alt_times)[win_alt]
                 model.set_dose(conc)
                 after = model.simulate(binding_params, times)[win]
-                model_milnes = np.append(model_milnes, after/before)
-                control = np.append(control, before)
+                model_milnes = np.append(model_milnes, after)
+                #control = np.append(control, before)
                 for i in range(swps*2-1):
                     if (alt_protocol is not None) & ((i % 2) == 0):
                         model.change_protocol(alt_protocol)
                         after = model.simulate(binding_params, alt_times, reset=False)[win_alt]
-                        model_milnes = np.append(model_milnes, after/before_alt)
-                        control	= np.append(control, before_alt)
+                        model_milnes = np.append(model_milnes, after)
+                        #control	= np.append(control, before_alt)
                         model.change_protocol(prot)
                     else:
                         after = model.simulate(binding_params, times, reset=False)[win]
-                        model_milnes = np.append(model_milnes, after/before)
-                        control = np.append(control, before)
+                        model_milnes = np.append(model_milnes, after)
+                        #control = np.append(control, before)
                 model_out[m][conc] = model_milnes
-                save_control = control
+                #save_control = control
             except:
                 model_out[m][conc] = np.ones(times.shape) * float('inf')
-                save_control = np.ones(times.shape) * float('inf')
-    return model_out, save_control, swps
+                #save_control = np.ones(times.shape) * float('inf')
+    return model_out, swps

protocols/Milnes_16102024_MA1_FP_RT.mmt

@@ -0,0 +1,6 @@
+[[protocol]]
+# Milnes protocol
+# Level Start           Length  Period  Multiplier
+-80     0.0             1350.0  25350   10
+0       1350.0          10000.0 25350   10
+-80     11350.0         14000.0 25350   10

src/fit_models.py

@@ -33,6 +33,7 @@ def parse_list_of_lists(s):
 parser.add_argument('-t', type=float, default=15e3, help='Max time')
 parser.add_argument('-b', type=parse_list_of_lists, default="[[1e3, 11e3]]", help='Protocol window(s) of interest')
 parser.add_argument('-c', type = str, help='Drug compound string')
+parser.add_argument('-d', action='store_true', help='Enable dual fitting of Milnes and optimal protocol data')
 args = parser.parse_args()
 
 def get_pars(model_num):
@@ -43,12 +44,21 @@ def get_pars(model_num):
         if herg_model != 'kemp' and herg_model != '2024_Joey_sis_25C':
             model = classes.ConcatMilnesModel(f'm{model_num}', protocol, times,
                                           win, conc, param_dict)
+            if args.d:
+                model_m = classes.ConcatMilnesModel(f'm{model_num}', 'protocols/Milnes_Phil_Trans.mmt', times_m,
+                                          win_m, conc, param_dict)
         elif herg_model == 'kemp':
             model = classes.ConcatMilnesModel(f'kemp-m{model_num}', protocol, times,
                                           win, conc, param_dict)
+            if args.d:
+                model_m = classes.ConcatMilnesModel(f'kemp-m{model_num}', 'protocols/Milnes_Phil_Trans.mmt', times_m,
+                                          win_m, conc, param_dict)
         elif herg_model == '2024_Joey_sis_25C':
             model = classes.ConcatMilnesModel(f'sis-m{model_num}', protocol, times,
                                           win, conc, param_dict)
+            if args.d:
+                model_m = classes.ConcatMilnesModel(f'sis-m{model_num}', 'protocols/Milnes_Phil_Trans.mmt', times_m,
+                                          win_m, conc, param_dict)
         # Load data
         u = np.loadtxt(
             f'{outdir}/fb_synthetic_conc_{conc}.csv',
@@ -57,10 +67,26 @@ def get_pars(model_num):
         )
         concat_time = u[:, 0]
         concat_milnes = u[:, 1]
+        if args.d:
+            # Load data
+            u_m = np.loadtxt(
+                f'{outdir.rsplit("/",1)[0]}/fb_synthetic_conc_{conc}.csv',
+                delimiter=',',
+                skiprows=1
+            )
+            concat_time_m = u_m[:, 0]
+            concat_milnes_m = u_m[:, 1]
         # Create single output problem
         problem = pints.SingleOutputProblem(model, concat_time, concat_milnes)
         likelihoods.append(classes.NormalRatioLogLikelihood(problem, mu_y))
-    f = pints.SumOfIndependentLogPDFs(likelihoods)
+        if args.d:
+            problem_m = pints.SingleOutputProblem(model_m, concat_time_m, concat_milnes_m)
+            likelihoods.append(classes.NormalRatioLogLikelihood(problem_m, mu_y_m))
+
+    if len(likelihoods) > 1:
+        f = pints.SumOfIndependentLogPDFs(likelihoods)
+    else:
+        f = likelihoods[0]
     bounds = boundaries.Boundaries(model_num, fix_hill=False, likelihood=True)
 
     # Fix random seed for reproducibility
@@ -114,8 +140,12 @@ def get_pars(model_num):
         concs = [30, 100, 300]
     elif args.c == 'quinidine':
         concs = [150, 500, 1500]
+    elif args.c == 'terfenadine':
+        concs = [30, 100, 300]
     elif args.c == 'verapamil':
         concs = [100, 300, 1000]
+    elif args.c == 'DMSO':
+        concs = [1]
     if herg_model != 'kemp' and herg_model != '2024_Joey_sis_25C':
         with open(f'methods/models/params/{p_path}', newline='') as csvfile:
             p_reader = csv.reader(csvfile)
@@ -133,9 +163,21 @@ def get_pars(model_num):
     win = np.zeros_like(conditions[0], dtype=bool)
     for condition in conditions:
         win |= condition
+    if args.d:
+        times_m = np.arange(0, 15e3, steps)
+        conditions_m = []
+        for b in [[1e3, 11e3]]:
+            conditions_m.append(((times_m >= b[0]) & (times_m < b[-1])))
+        win_m = np.zeros_like(conditions_m[0], dtype=bool)
+        for condition in conditions_m:
+            win_m |= condition
     # read fitted splines
     dfy = pd.read_csv(f"{outdir}/synth_Y_fit.csv")
     mu_y = np.array(dfy['0'])
+    if args.d:
+        # read fitted splines
+        dfy_m = pd.read_csv(f"{outdir.rsplit('/',1)[0]}/synth_Y_fit.csv")
+        mu_y_m = np.array(dfy_m['0'])
     # fit model
     pars, sc = get_pars(model_num)
     print(f'{model_num}: {pars}, {sc}')

src/fit_spline.py

@@ -30,6 +30,9 @@ def main(input, output, lambda_, t_steps, data, max_t):
         t_swp = t_steps[j]
         t_total = 0
         swps = int(np.floor(250000/max_t))*6
+    elif data == 'milnes real':
+        t_swp = 10000
+        swps = 20
     else:
         lens = [3340, 3330, 3340, 3330, 3330]
         t_total = 0
@@ -39,7 +42,7 @@ def main(input, output, lambda_, t_steps, data, max_t):
 
     # Fit splines
     for i in range(swps):
-        if data == 'milnes':
+        if data == 'milnes' or data == 'milnes real':
             df_rep = df_all[(df_all['t'] >= t_swp * i) & (df_all['t'] < t_swp * (i + 1))][['t', 'x']]
             knots = np.arange(t_swp * i, t_swp * (i + 1) + t_swp/2, t_swp/2)
         elif data == 'opt':

src/optimise_protocol.py

@@ -22,6 +22,7 @@
 parser.add_argument('-e', type=str, default='joey_sis', help='hERG model parameters')
 parser.add_argument('-o', type=str, required=True, help='output folder for synthetic data')
 parser.add_argument('-c', type=str, help='drug compound')
+parser.add_argument('-r', action='store_true', help='set true for real data')
 args = parser.parse_args()
 
 def get_opt_prot(model_pars, herg, v_steps, t_steps, p0, CMAES_pop = 10, max_iter = 5, alt_protocol = None):
@@ -49,32 +50,32 @@ def __call__(self, p):
                 prot = funcs.create_protocol(self.v_st, self.t_st)
             # generate model output under the new protocol
             if alt_protocol is not None:
-                model_out, cont, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs,
+                model_out, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs,
                                           alt_protocol = prot_alt, alt_times = np.arange(0, np.sum(self.t_st_alt), 10), wins = [1e3, np.sum(self.t_st[1:4])],
                                           wins_alt = [1e3, np.sum(self.t_st_alt[1:4])])
             else:
-                model_out, cont, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs, wins = [1e3, np.sum(self.t_st[1:4])])
+                model_out, swps = funcs.model_outputs(model_pars, herg, prot, times = np.arange(0, np.sum(self.t_st), 10), concs = concs, wins = [1e3, np.sum(self.t_st[1:4])])
             ### loop through models and concentrations to get traces
             all_traces = []
             for m in model_out:
                 model_trace = []
                 for c in concs:
                     model_trace = np.append(model_trace, model_out[m][c])
                 all_traces.append(model_trace)
-            conts = []
-            for c in concs:
-                conts = np.append(conts, cont)
-            #ssq = [sum((m-n)**2) for i,m in enumerate(all_traces) for j,n in enumerate(all_traces) if i < j]
+            #conts = []
+            #for c in concs:
+            #    conts = np.append(conts, cont)
+            ssq = [sum((m-n)**2) for i,m in enumerate(all_traces) for j,n in enumerate(all_traces) if i < j]
             # calculate expected likelihood ratio (assuming normal ratio data)
-            lhoods = []
-            for i,m in enumerate(all_traces):
-                for j,n in enumerate(all_traces):
-                    if i < j:
-                        top = np.exp(-conts**2*((m**2+1)/(2*sd**2))) + np.sqrt(np.pi/2)*(conts/sd)*np.sqrt(1+m**2)*special.erf((conts/(np.sqrt(2)*sd))*np.sqrt(1+m**2))
-                        bottom = np.exp(-conts**2*((n**2+1)/(2*sd**2))) + np.sqrt(np.pi/2)*(conts/sd)*((1+m*n)/np.sqrt(1+m**2))*special.erf((conts*(1+m*n))/(np.sqrt(2)*sd*np.sqrt(1+m**2)))*np.exp(-(conts**2*(m-n)**2)/(2*sd**2*(1+m**2)))
-                        lhoods.append(-2*np.sum(np.log(top/bottom)))
-            #out = -np.median(ssq)
-            out = np.median(lhoods)
+            #lhoods = []
+            #for i,m in enumerate(all_traces):
+            #    for j,n in enumerate(all_traces):
+            #        if i < j:
+            #            top = np.exp(-conts**2*((m**2+1)/(2*sd**2))) + np.sqrt(np.pi/2)*(conts/sd)*np.sqrt(1+m**2)*special.erf((conts/(np.sqrt(2)*sd))*np.sqrt(1+m**2))
+            #            bottom = np.exp(-conts**2*((n**2+1)/(2*sd**2))) + np.sqrt(np.pi/2)*(conts/sd)*((1+m*n)/np.sqrt(1+m**2))*special.erf((conts*(1+m*n))/(np.sqrt(2)*sd*np.sqrt(1+m**2)))*np.exp(-(conts**2*(m-n)**2)/(2*sd**2*(1+m**2)))
+            #            lhoods.append(-2*np.sum(np.log(top/bottom)))
+            out = -np.median(ssq)
+            #out = np.median(lhoods)
             return out
 
         def n_parameters(self):
@@ -99,14 +100,14 @@ def n_parameters(self):
     transformation = pints.RectangularBoundariesTransformation(boundaries)
 
     # define initial standard deviation around voltages and times during optimisation
-    step_v = [5]*v_steps.count(np.nan)
-    step_t = [50]*t_steps.count(np.nan)
+    step_v = [10]*v_steps.count(np.nan)
+    step_t = [500]*t_steps.count(np.nan)
 
     if alt_protocol is not None:
         design = DrugBind(p0+alt_protocol)
     else:
         design = DrugBind(p0)
-    
+
     # Fix random seed for reproducibility
     np.random.seed(101)
 
@@ -117,12 +118,14 @@ def n_parameters(self):
         q0 = boundaries.sample()[0]
         try:
             temp = design(q0)
+            print(f'temp score: {temp}')
             if temp < score:
                 score = temp
                 q0save = q0
         except Exception as e:
             print(f"An error occurred: {e}")
 
+    print(f'init_score: {score}')
     # Define optimiser
     optimiser = pints.CMAES
     if alt_protocol is not None:
@@ -211,6 +214,14 @@ def main(model_nums, max_time, bounds, herg, output_folder):
         writer.writerow(t_steps)
 
 if __name__ == "__main__":
-    concs = parameters.drug_concs[args.c]
+    if args.r:
+        if args.c == 'verapamil':
+            concs=[100,300,1000]
+        elif args.c == 'bepridil':
+            concs=[30,100,300]
+        elif args.c == 'terfenadine':
+            concs=[100,300,1000]
+    else:
+        concs = parameters.drug_concs[args.c]
     m_list = ast.literal_eval(args.m)
     main(m_list, args.t, args.b, args.e, args.o)