diff --git a/Fig3_S3/Growth_inh_density.py b/Fig3_S3/Growth_inh_density.py
deleted file mode 100644
index d945c1b..0000000
--- a/Fig3_S3/Growth_inh_density.py
+++ /dev/null
@@ -1,190 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Oct 19 11:07:03 2023
-
-@author: nicagutu
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import pandas as pd
-from scipy.signal import savgol_filter
-from scipy.optimize import curve_fit
-from pyboat import WAnalyzer
-from pyboat import ensemble_measures as em
-import os
-import seaborn as sns
-    
-plt.rcParams.update({'font.size': 28})
-plt.rcParams['svg.fonttype'] = 'none'    
-  
-def exponential_func(t, N0, r):
-    K = max(ydata_smoothed)
-    return K/(1+((K-N0)/N0)*np.exp(-r*t))
-
-def linear_func(t, a, b):
-    return a*t+b
-
-def growthcurve(x0, t, k, kl, LC50, Sm, c, SC50, h): #cell count
-    return x0*np.exp(t*k*(1-(Sm*c**h)/(SC50**h+c**h))-t*(kl*c**h)/(LC50**h+c**h))
-
-def exp_func(b, x):
-    return 0.068*np.exp(b*x)
-
-def exp_func2(b, x):
-    return 0.13*np.exp(-b*x)
-
-dt = 1 
-lowT = 16
-highT = 40
-periods = np.linspace(lowT, highT, 200)
-wAn = WAnalyzer(periods, dt, time_unit_label='hours')
-            
-#%%Data
-path = 'RawData/'
-output = 'Figures/'
-
-file = 'Circadian_growth_TGFbeta.xlsx'
-
-data1_confl = pd.read_excel(path+file, sheet_name=0, skiprows=[0], header=[0])
-data1_green = pd.read_excel(path+file, sheet_name=1, skiprows=[0], header=[0]) 
-
-data2_confl = pd.read_excel(path+file, sheet_name=4, skiprows=[0], header=[0])
-data2_green = pd.read_excel(path+file, sheet_name=5, skiprows=[0], header=[0]) 
-
-#%%Average results for confluence
-data1_confl = data1_confl.drop(columns=['Date Time', 'Elapsed'])
-data2_confl = data2_confl.drop(columns=['Date Time', 'Elapsed'])
-data_confl = (data1_confl+data2_confl)/2
-data_confl = data_confl.drop([2])
-
-line_colors = ['tab:green', 'tab:orange', 'tab:blue']
-labels1 = ['10uM', '5uM', 'untreated']
-labels2 = ['high', 'medium', 'low']
-
-columns1 = ['A1', 'A2', 'A6']
-columns2 = ['A6', 'B6', 'C6']
-
-fit_params1 = []
-cov_error1 = []
-
-doubling_times1 = []
-
-plt.figure(figsize=(12,10))
-for col, label, color in zip(columns1, labels1, line_colors):
-    minimum = data_confl[col].min()
-    index = data_confl[col].idxmin()
-    ydata = data_confl[col][index:]#/minimum
-    
-    #shorter range of data for linear fit of the first 24h
-    time = ydata.index.values#[0:24]
-    ydata = ydata#[0:24]
-
-    ydata_smoothed = savgol_filter(ydata, window_length=9, polyorder=2)
-    params, covariance = curve_fit(exponential_func, time, ydata_smoothed)
-    # params, covariance = curve_fit(linear_func, time, ydata_smoothed)
-
-    errors = np.sqrt(np.diag(covariance))
-    y_fit = exponential_func(time, *params)
-    # y_fit = linear_func(time, *params)
-    
-    var = 1
-    fit_params1.append(params[var])
-    cov_error1.append(errors[var])
-    doubling_times1.append(np.log(2)/params[var])
-    
-    plt.plot(time, ydata_smoothed, label=label, linewidth=5, color=color)
-    plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, color=color, alpha=0.5)
-    # plt.plot(time, y_fit, '--', label='fit: {:.2f} $\pm$ {:.2f}'.format(params[0], errors[0]), linewidth=5, color=color, alpha=0.5)
-
-plt.ylabel('Confluence (%)')
-plt.xlabel('Time [hours]')
-# plt.xticks([0,24,48,72,96,120])
-plt.legend(loc='best')
-# plt.savefig(output+'Confluence_change_density.svg')
-plt.show()
-
-fit_params2 = []
-cov_error2 = []
-
-doubling_times2 = []
-
-plt.figure(figsize=(12,10))
-for col, label, color in zip(columns2, labels2, line_colors):
-    minimum = data_confl[col].min()
-    index = data_confl[col].idxmin()
-    ydata = data_confl[col][index:]#/minimum
-    
-    #shorter range of data for linear fit of the first 24h
-    time = ydata.index.values#[0:24]
-    ydata = ydata#[0:24]
-
-    ydata_smoothed = savgol_filter(ydata, window_length=9, polyorder=2)
-    params, covariance = curve_fit(exponential_func, time, ydata_smoothed)
-    # params, covariance = curve_fit(linear_func, time, ydata_smoothed)
-
-    errors = np.sqrt(np.diag(covariance))
-    y_fit = exponential_func(time, *params)
-    # y_fit = linear_func(time, *params)
-
-    var = 1
-    fit_params2.append(params[var])
-    cov_error2.append(errors[var])
-    doubling_times2.append(np.log(2)/params[var])
-    
-    plt.plot(time, ydata_smoothed, label=label, linewidth=5, color=color)
-    plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, color=color, alpha=0.5)
-    # plt.plot(time, y_fit, '--', label='fit: {:.2f} $\pm$ {:.2f}'.format(params[0], errors[0]), linewidth=5, color=color, alpha=0.5)
-
-plt.ylabel('Confluence (%)')
-plt.xlabel('Time [hours]')
-# plt.xticks([0,24,48,72,96,120])
-plt.legend(loc='best')
-# plt.savefig(output+'Confluence_change_density.svg')
-plt.show()
-
-fit_params1 = np.array(fit_params1)
-cov_error1 = np.array(cov_error1)
-print(fit_params1)
-print(cov_error1)
-
-p0 = fit_params1[-1]
-normalized_params = fit_params1 / p0
-normalized_errors = normalized_params * np.sqrt((cov_error1 / fit_params1) ** 2 + (cov_error1[-1] / p0) ** 2)
-
-plt.figure(figsize=(10,8))
-plt.errorbar(labels1, normalized_params, yerr=normalized_errors, fmt='o')
-plt.xlabel('Inhibitor concentration') #Inhibitor concentration
-plt.ylabel('Relative growth rate')
-plt.ylim([0.5, 1.1])
-plt.savefig(output+'Logistic_fit_growth_curves_inhibitor.svg')
-# plt.savefig(output+'Linear_fit_24h_growth_curves_inhibitor.svg')
-plt.show()
-
-fit_params2 = np.array(fit_params2)
-cov_error2 = np.array(cov_error2)
-print(fit_params2)
-print(cov_error2)
-
-p0 = fit_params2[0]
-normalized_params = fit_params2 / p0
-normalized_errors = normalized_params * np.sqrt((cov_error2 / fit_params2) ** 2 + (cov_error2[0] / p0) ** 2)
-
-plt.figure(figsize=(10,8))
-plt.errorbar(labels2, normalized_params, yerr=normalized_errors, fmt='o')
-plt.xlabel('Seeding density') #Inhibitor concentration
-plt.ylabel('Relative growth rate')
-plt.ylim([0.5, 1.1])
-plt.savefig(output+'Logistic_fit_growth_curves_density.svg')
-# plt.savefig(output+'Linear_fit_24h_growth_curves_density.svg')
-plt.show()
-
-
-    
-    
-    
-    
-
-
-
diff --git a/Fig3_S3/Norm_relative_growth_WT_dKO.py b/Fig3_S3/Norm_relative_growth_WT_dKO.py
deleted file mode 100644
index fd133fc..0000000
--- a/Fig3_S3/Norm_relative_growth_WT_dKO.py
+++ /dev/null
@@ -1,139 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Feb 27 13:53:06 2024
-
-@author: nicagutu
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import pandas as pd
-from scipy.signal import savgol_filter
-from scipy.optimize import curve_fit
-from pyboat import WAnalyzer
-    
-plt.rcParams.update({'font.size': 24})
-plt.rcParams['svg.fonttype'] = 'none'      
-
-def exponential_func(t, N0, r):
-    K = max(ydata)
-    return K/(1+((K-N0)/N0)*np.exp(-r*t))
-
-
-dt = 1.5 # the sampling interval, 0.5hours
-lowT = 16
-highT = 32
-periods = np.linspace(lowT, highT, 200)
-wAn = WAnalyzer(periods, dt, time_unit_label='hours')
-
-path = 'Raw_data/'
-output = '...'
-file1 = 'U2OS_3osc_dko_1stplate.xlsx'
-file2 = 'U2OS_3osc_dko_3rdplate.xlsx'
-
-data1_confl = pd.read_excel(path+file1, sheet_name=0, skiprows=[0], header=[0])
-data1_green = pd.read_excel(path+file1, sheet_name=1, skiprows=[0], header=[0]) 
-
-data2_confl = pd.read_excel(path+file2, sheet_name=0, skiprows=[0], header=[0])
-data2_green = pd.read_excel(path+file2, sheet_name=1, skiprows=[0], header=[0]) 
-
-#%%Average confluence results 1st plate
-data1_confl = data1_confl.drop(columns=['Date Time', 'Elapsed'])
-
-index = (max(data1_confl.idxmin()))
-
-data1_confl_smoothed = pd.DataFrame()
-
-for col in data1_confl:
-    ydata = data1_confl[col][index:]
-    data1_confl_smoothed[col] = savgol_filter(ydata/ydata[index], window_length=9, polyorder=2)
-
-time = data1_confl_smoothed.index.values*1.5+index
-reset = 57#-index*1.5
-    
-aver_dKo_high_dmso = data1_confl_smoothed[['A1','A2','A3']].mean(axis=1) 
-aver_dKo_high_inh = data1_confl_smoothed[['A4', 'A5', 'A6']].mean(axis=1) 
-
-std_dKo_high_dmso = data1_confl_smoothed[['A1', 'A2', 'A3']].std(axis=1) 
-std_dKo_high_inh = data1_confl_smoothed[['A4', 'A5', 'A6']].std(axis=1) 
-
-aver_3osc_high_dmso = data1_confl_smoothed[['C1', 'C2', 'C3']].mean(axis=1)
-aver_3osc_high_inh = data1_confl_smoothed[['C4', 'C5', 'C6']].mean(axis=1) 
-
-std_3osc_high_dmso = data1_confl_smoothed[['C1', 'C2','C3']].std(axis=1) 
-std_3osc_high_inh = data1_confl_smoothed[['C4', 'C5','C6']].std(axis=1) 
-
-plt.figure(figsize=(10,8))
-
-plt.plot(time, aver_dKo_high_dmso, linewidth=5, color='blue', label='dKO high DMSO')
-plt.errorbar(time, aver_dKo_high_dmso, yerr=std_dKo_high_dmso, color='blue', alpha=0.3)
-# ydata = aver_dKo_high_dmso
-# params, covariance = curve_fit(exponential_func, time, ydata)
-# y_fit = exponential_func(time, *params)
-# print('r', params[1])
-# errors = np.sqrt(np.diag(covariance))
-# print('error', errors[1])
-# plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-
-plt.plot(time, aver_dKo_high_inh, '--', linewidth=5, color='tab:green', label='dKO high inh')
-plt.errorbar(time, aver_dKo_high_inh, yerr=std_dKo_high_inh, color='tab:green', alpha=0.3)
-# ydata = aver_dKo_high_inh
-# params, covariance = curve_fit(exponential_func, time, aver_dKo_high_inh)
-# y_fit = exponential_func(time, *params)
-# print('r', params[1])
-# errors = np.sqrt(np.diag(covariance))
-# print('error', errors[1])
-# plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-
-plt.plot(time, aver_3osc_high_dmso, linewidth=5, color='brown', label='wt high DMSO')
-plt.errorbar(time, aver_3osc_high_dmso, yerr=std_3osc_high_dmso, color='brown', alpha=0.3)
-# ydata = aver_3osc_high_dmso
-# params, covariance = curve_fit(exponential_func, time, ydata)
-# y_fit = exponential_func(time, *params)
-# print('r', params[1])
-# errors = np.sqrt(np.diag(covariance))
-# print('error', errors[1])
-# plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-
-plt.plot(time, aver_3osc_high_inh, '--', linewidth=5, color='tab:orange', label='wt high inh')
-plt.errorbar(time, aver_3osc_high_inh, yerr=std_3osc_high_inh, color='tab:orange', alpha=0.3)
-# ydata = aver_3osc_high_inh
-# params, covariance = curve_fit(exponential_func, time, ydata)
-# y_fit = exponential_func(time, *params)
-# print('r', params[1])
-# errors = np.sqrt(np.diag(covariance))
-# print('error', errors[1])
-# plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-
-plt.xlabel('Time [hours]')
-plt.ylabel('Normalized confluence')
-plt.legend(loc='best')
-plt.xticks([24,48,72,96,120,148])
-plt.savefig(output+'Raw_confl_3osc_dko_DMSO_inh_high_dens.svg')
-plt.show()
-
-
-plt.figure(figsize=(10,8))
-
-plt.plot(time, aver_dKo_high_dmso/aver_dKo_high_inh, linewidth=5, label='dKO high')
-ydata = aver_dKo_high_dmso/aver_dKo_high_inh
-params, covariance = curve_fit(exponential_func, time, ydata)
-y_fit = exponential_func(time, *params)
-plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-errors = np.sqrt(np.diag(covariance))
-print(errors[1])
-
-plt.plot(time, aver_3osc_high_dmso/aver_3osc_high_inh, linewidth=5, label='3osc high')
-ydata = aver_3osc_high_dmso/aver_3osc_high_inh
-params, covariance = curve_fit(exponential_func, time, ydata)
-y_fit = exponential_func(time, *params)
-plt.plot(time, y_fit, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5)
-errors = np.sqrt(np.diag(covariance))
-print(errors[1])
-
-plt.xlabel('Time [hours]')
-plt.ylabel('Relative confluence')
-plt.legend(loc='best')
-plt.savefig(output+'Relative_confl_3osc_dko_DMSO_inh_high_dens_reset0.svg')
-plt.show()
diff --git a/Fig3_S3/cell_growth_plot.py b/Fig3_S3/cell_growth_plot.py
deleted file mode 100644
index 9b65b8e..0000000
--- a/Fig3_S3/cell_growth_plot.py
+++ /dev/null
@@ -1,119 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 30 11:52:35 2021
-
-@author: gutunn
-"""
-
-import pandas as pd
-import matplotlib.pyplot as plt
-import numpy as np
-from scipy.signal import savgol_filter
-from scipy.optimize import curve_fit
-
-plt.rcParams.update({'font.size':24})
-plt.rcParams['svg.fonttype'] = 'none'
-
-def exponential_func(t, N0, r):
-    K = ydata_smoothed[-1]
-    return K/(1+((K-N0)/N0)*np.exp(-r*t))
-
-def linear_func(t, a, b):
-    return a*t+b
-
-
-table_path = 'RawData'
-output = 'Figures/'
-
-
-condition_combine_dict = {101:'high_density_untreated', 21:'high_density_5uM',1:'high_density_10uM'}
-# condition_combine_dict = {101:'high_density_untreated', 141:'medium_density_untreated',341:'low_density_untreated'}
-
-# count cell numbers for each condition over frames, and save into a dictionary
-count_dict = {}
-for j in condition_combine_dict:    
-    count_list = [0 for x in range(232)]
-    for i in range(20):
-        df = pd.read_csv('%s/xy%03d/xy%03d-t_tracking_table.csv' % (table_path, j+i, j+i))
-        
-        count_list_pos = []
-        for k in range(232):
-            c = len(set(df.loc[df['frame']==k]['trackId']))
-            count_list_pos.append(c)
-            
-        count_list = [(l+m) for l,m in zip(count_list,count_list_pos)]
-    
-    count_list = [x/count_list[0] for x in count_list] # this line defines how you normalize cell counts
-    count_dict[j] = count_list
-
-
-# plot results
-density = 1
-dose = 3
-
-pos_list = [ x for x in condition_combine_dict]
-color_list = ['gold','tab:green','tab:blue']
-label_list = ['untreated','5uM','10uM'] #'0uM','1.25uM','2.5uM',
-# label_list = ['high density','50% density','25% density']
-
-print(pos_list)
-
-fit_params1 = []
-cov_error1 = []
-doubling_times1 = []
-
-fig, ax = plt.subplots()
-for j in range(len(label_list)):
-    index = pos_list[j]
-    ydata = count_dict[index]
-
-    time = np.arange(0,116,0.5)   
-    ydata_smoothed = savgol_filter(ydata, window_length=5, polyorder=3)
-
-    tf = np.argmax(ydata_smoothed[0:150])
-    time_fit = time[10:tf]
-    ydata_fit = ydata_smoothed[10:tf]
-    # params, covariance = curve_fit(exponential_func, time_fit, ydata_fit)
-    params, covariance = curve_fit(linear_func, time_fit, ydata_fit)
-
-    errors = np.sqrt(np.diag(covariance))
-    # y_fit = exponential_func(time, *params)
-    y_fit = linear_func(time_fit, *params)
-    
-    ax.plot(time, ydata_smoothed, label = label_list[j], color = color_list[j],linewidth=5)
-    ax.plot(time_fit, y_fit, '--', label='fit: {:.4f} $\pm$ {:.4f}'.format(params[0], errors[0]), alpha=0.5, color = color_list[j],linewidth=5)
-
-    var = 0 # 0 for linear fit and 1 for exp fit
-    fit_params1.append(params[var])
-    cov_error1.append(errors[var])
-    doubling_times1.append(np.log(2)/params[var])
-
-ax.axvline(x=74, color='grey', linestyle='--', alpha=0.5)
-ax.legend()
-ax.set_xlabel('Time [hours]')
-ax.set_ylabel('Normalized number of cells')
-fig.set_figheight(12)
-fig.set_figwidth(14)
-# plt.savefig(output+'Cell_counts_changing_density.svg')
-plt.show()
-
-fit_params1 = np.array(fit_params1)
-cov_error1 = np.array(cov_error1)
-print(fit_params1)
-print(cov_error1)
-
-p0 = fit_params1[0]
-normalized_params = fit_params1 / p0
-normalized_errors = normalized_params * np.sqrt((cov_error1 / fit_params1) ** 2 + (cov_error1[0] / p0) ** 2)
-
-plt.figure(figsize=(10,8))
-plt.errorbar(label_list, normalized_params, yerr=normalized_errors, fmt='o')
-plt.xlabel('Inhibitor concentration') #Inhibitor concentration
-plt.ylabel('Growth rate')
-# plt.ylim([0.5, 1.1])
-# plt.savefig(output+'Logistic_fit_growth_curves_inhibitor.svg')
-# plt.savefig(output+'Linear_fit_until_saturation_obj_count_density.svg')
-plt.show()
-
-print(doubling_times1) 
-    
\ No newline at end of file
diff --git a/Fig3_S3/imt_cond_total_num_divisions.py b/Fig3_S3/imt_cond_total_num_divisions.py
deleted file mode 100644
index 1177469..0000000
--- a/Fig3_S3/imt_cond_total_num_divisions.py
+++ /dev/null
@@ -1,124 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Oct 11 12:50:35 2024
-
-@author: nicagutu
-"""
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Sep 24 12:55:32 2024
-
-@author: nicagutu
-"""
-
-import pandas as pd
-import numpy as np
-import matplotlib.pyplot as plt
-from pyboat import WAnalyzer
-from pyboat import ensemble_measures as em
-from pyboat import plotting as pl 
-import seaborn as sns
-import random
-from scipy.stats import kurtosis
-import scipy.stats as stats
-from scipy.stats import ttest_ind
-
-plt.rcParams.update({'font.size': 22})    
-plt.rcParams['svg.fonttype'] = 'none'
-
-path = 'RawData/'
-output = 'Figures/'
-
-#%% Changing TGF-beta inhibitor dose
-dose = ['untreated', '5uM', '10uM']
-density = 'high'#['high', 'medium', 'low']
-
-# Create an empty dictionary to store intermitotic times (IMT) for each division group across all doses
-imt_all_conditions = {str(num): {dose_level: [] for dose_level in dose} for num in range(2, 6)}
-
-for i in dose:
-    condition = f'{density}_density_{i}_'
-    print(condition)
-
-    # Load data
-    data = pd.read_excel(path + f'divisions_{condition}.xlsx')
-
-    num_divisions = []
-    for col in data:
-        div_events = data[col]
-        division_times = np.where(div_events == 1)[0]
-        
-        num_divisions.append(len(division_times))
-        if 2 <= len(division_times) <= 5:
-            imt_single = []
-            for ii in range(len(division_times) - 1):
-                imt_single.append((division_times[ii + 1] - division_times[ii]) * 0.5)
-            # Append the average IMT for this particular cell to the corresponding division group for this dose
-            imt_all_conditions[str(len(division_times))][i].append(np.mean(imt_single))
-
-# Now let's plot boxplots per division group, with different conditions in each boxplot
-fig, axes = plt.subplots(2, 2, figsize=(14, 12))
-axes = axes.flatten()
-
-for idx, num_div in enumerate(range(2, 6)):
-    # Collect the data for the current division group
-    division_group_data = [imt_all_conditions[str(num_div)][dose_level] for dose_level in dose]
-    
-    # Create a boxplot
-    axes[idx].boxplot(division_group_data, labels=dose, showfliers=False)
-    axes[idx].set_title(f'{num_div} divisions')
-    axes[idx].set_ylim([15, 50])
-    axes[idx].set_ylabel('IMT [hours]')
-
-plt.tight_layout()
-plt.savefig(output+'IMT_dose_condition_total_divisions.svg')
-plt.show()    
-
-#%% Changing seeding density
-
-dose = 'untreated'#['untreated', '5uM', '10uM']
-density = ['high', 'medium', 'low']
-
-# Create an empty dictionary to store intermitotic times (IMT) for each division group across all doses
-imt_all_conditions = {str(num): {dens_level: [] for dens_level in density} for num in range(2, 6)}
-
-for i in density:
-    condition = f'{i}_density_{dose}_'
-    print(condition)
-
-    # Load data
-    data = pd.read_excel(path + f'divisions_{condition}.xlsx')
-
-    num_divisions = []
-    for col in data:
-        div_events = data[col]
-        division_times = np.where(div_events == 1)[0]
-        
-        num_divisions.append(len(division_times))
-        if 2 <= len(division_times) <= 5:
-            imt_single = []
-            for ii in range(len(division_times) - 1):
-                imt_single.append((division_times[ii + 1] - division_times[ii]) * 0.5)
-            # Append the average IMT for this particular cell to the corresponding division group for this dose
-            imt_all_conditions[str(len(division_times))][i].append(np.mean(imt_single))
-
-# Now let's plot boxplots per division group, with different conditions in each boxplot
-fig, axes = plt.subplots(2, 2, figsize=(14, 12))
-axes = axes.flatten()
-
-for idx, num_div in enumerate(range(2, 6)):
-    # Collect the data for the current division group
-    division_group_data = [imt_all_conditions[str(num_div)][dens_level] for dens_level in density]
-    
-    # Create a boxplot
-    axes[idx].boxplot(division_group_data, labels=density, showfliers=False)
-    axes[idx].set_title(f'{num_div} divisions')
-    axes[idx].set_ylim([15, 50])
-    axes[idx].set_ylabel('IMT [hours]')
-
-plt.tight_layout()
-plt.savefig(output+'IMT_density_condition_total_divisions.svg')
-plt.show()