diff --git a/Fig4_S4/CC_model/CC_illu.py b/Fig4_S4/CC_model/CC_illu.py deleted file mode 100644 index 4101936..0000000 --- a/Fig4_S4/CC_model/CC_illu.py +++ /dev/null @@ -1,157 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -# plt.style.use('K_PAPER') -from tools import * -from scipy.signal import find_peaks -import tqdm -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from poincare_oscillator import coupling -from multiprocessing import Pool -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' - -plt.rcParams.update({'font.size': 14}) -plt.rcParams['svg.fonttype'] = 'none' - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def log2(x, xhalf, beta): - return 1/(1 + np.exp(beta*(x - xhalf))) - -def run(circ, period): - # circ = np.roll(circ, np.random.randint(30,200,1)[0]) - RT = 0 - circ = (1/2 + 1/2*circ)**2 - # circ /= np.max(circ) - circ = np.nan_to_num(circ) - - Vol = 2 - - G1 = 3 - G2 = 0 - S = 0 - Div = 0 - AM = 0 - # RT = 0 - Sav = np.zeros((100000, 6)) - - k1 = 50 - k2 = 20 - k3 = 200 - k4 = 9 - k5 = 0.7 - # k6 = 0.5 - k6 = 0.9 - # k6 = 0.01 - k7 = 1 - k8 = 5e-3 - - alpha = .9 - # alpha = 0 - - GMAX = period*3.8 + stats.expon.rvs(0, 1/.007,1)[0] - GMAX = int(GMAX) - - - - rSav = np.zeros((4,280)) - - for i in range(100000): - - ptime = int((RT/.1)) - - dG1p = k1*(circ[ptime]*alpha + (1-alpha))*(1-S)*Vol - dG1n = k2 * G1 / (G1 + k3)*Vol - dS = k5*max(G1-GMAX*Vol,0)*(1-S)#*circ[ptime] - - dG2p = k4*S*Vol - dDiv = k6 * log2(G2, 150*Vol, -1)*Vol#*circ[ptime] - dAM = (1-S)*k7*Vol - dA = 0*k8*log2(AM, 10*Vol, -1)*Vol - - rates = np.array([dG1p, dG1n, dG2p, dDiv, dS, dAM, dA]) - trates = np.sum(rates) - DT = 1/trates * np.log(1/np.random.uniform(0,1,1)[0]) - rates /= trates - - ran = np.random.uniform(0,1,1)[0] - tsum = 0 - idx = 0 - for j in range(7): - tsum += rates[j] - if ran < tsum: - idx = j - break - - if idx == 0: - G1 += 1 - - elif idx == 1: - G1 -= 1 - - elif idx == 2: - G2 += 1 - - elif idx == 3: - Div += 1 - G2 = 0 - S = 0 - rSav[0,int(RT*2):] = Div - - elif idx == 4: - S = 1 - AM = 1 - - elif idx == 5: - AM += 1 - - elif idx == 6: - Sav[i:, :] = np.nan - rSav[-1,int(RT*2)] = 1 - - break - - Sav[i,0] = G1 - Sav[i,1] = G2 - Sav[i,2] = circ[ptime] - Sav[i,3] = RT - Sav[i,4] = dA - Sav[i,5] = AM - rSav[1,int(RT*2)] = G1 - rSav[2,int(RT*2)] = G2 - - - RT += DT - - if RT > 140: - Sav[i:, :] = np.nan - break - fig, ax = plt.subplots(1,1, figsize = (6,4)) - ax.plot(Sav[000:3500,3], Sav[1000:4500,0]) - ax.plot(Sav[000:3500,3], Sav[1000:4500,1]) - ax.set(ylabel = 'Phase intensity [AU]', xlabel = 'Time[h]') - # ax.plot(Sav[:,1]) - ax.legend(['G1', 'S/G2/M'], title = 'Phase') - ax.set_xticks(np.arange(0,84,24), np.arange(0,84,24)) - plt.tight_layout() - fig.savefig('sim_trace.svg') - plt.show() - - return Sav, rSav - - -i = 100 -data = np.genfromtxt(f'nica_output/circ_none.csv', delimiter = ',') -period_cc = np.genfromtxt(f'nica_output/circ_dist_none.csv', delimiter = ',') -Sav, rSav = run(data[i,:], period_cc[i]) - - - - - - - - diff --git a/Fig4_S4/CC_model/G1_coupled_system.py b/Fig4_S4/CC_model/G1_coupled_system.py deleted file mode 100644 index a2d90aa..0000000 --- a/Fig4_S4/CC_model/G1_coupled_system.py +++ /dev/null @@ -1,207 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('K_PAPER') -from tools import * -from scipy.signal import find_peaks -import tqdm -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from poincare_oscillator import coupling -from multiprocessing import Pool -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' - - -''' -author @malthenielse -malthe.nielsen@nbi.ku.dk -''' - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def log2(x, xhalf, beta): - return 1/(1 + np.exp(beta*(x - xhalf))) - -def run(circ, period): - # circ = np.roll(circ, np.random.randint(30,200,1)[0]) - RT = 0 - circ = (1/2 + 1/2*circ)**2 - # circ /= np.max(circ) - circ = np.nan_to_num(circ) - - Vol = 1 - - G1 = 3 - G2 = 0 - S = 0 - Div = 0 - AM = 0 - # RT = 0 - Sav = np.zeros((100000, 6)) - - k1 = 50 - k2 = 20 - k3 = 200 - k4 = 9 - k5 = 0.7 - # k6 = 0.5 - k6 = 0.9 - # k6 = 0.01 - k7 = 1 - k8 = 5e-3 - - alpha = .9 - # alpha = 0 - - GMAX = period*3.8 + stats.expon.rvs(0, 1/.007,1)[0] - GMAX = int(GMAX) - - - - rSav = np.zeros((4,280)) - - for i in range(100000): - - ptime = int((RT/.1)) - - dG1p = k1*(circ[ptime]*alpha + (1-alpha))*(1-S)*Vol - dG1n = k2 * G1 / (G1 + k3)*Vol - dS = k5*max(G1-GMAX*Vol,0)*(1-S)#*circ[ptime] - - dG2p = k4*S*Vol - dDiv = k6 * log2(G2, 150*Vol, -1)*Vol#*circ[ptime] - dAM = (1-S)*k7*Vol - dA = 0*k8*log2(AM, 10*Vol, -1)*Vol - - rates = np.array([dG1p, dG1n, dG2p, dDiv, dS, dAM, dA]) - trates = np.sum(rates) - DT = 1/trates * np.log(1/np.random.uniform(0,1,1)[0]) - rates /= trates - - ran = np.random.uniform(0,1,1)[0] - tsum = 0 - idx = 0 - for j in range(7): - tsum += rates[j] - if ran < tsum: - idx = j - break - - if idx == 0: - G1 += 1 - - elif idx == 1: - G1 -= 1 - - elif idx == 2: - G2 += 1 - - elif idx == 3: - Div += 1 - G2 = 0 - S = 0 - rSav[0,int(RT*2):] = Div - - elif idx == 4: - S = 1 - AM = 1 - - elif idx == 5: - AM += 1 - - elif idx == 6: - Sav[i:, :] = np.nan - rSav[-1,int(RT*2)] = 1 - - break - - Sav[i,0] = G1 - Sav[i,1] = G2 - Sav[i,2] = circ[ptime] - Sav[i,3] = RT - Sav[i,4] = dA - Sav[i,5] = AM - rSav[1,int(RT*2)] = G1 - rSav[2,int(RT*2)] = G2 - - - RT += DT - - if RT > 140: - Sav[i:, :] = np.nan - break - - return Sav, rSav - -# coupling(0.0000, 'none') -# coupling(0.0005, 'medium') -# coupling(0.001, 'high') -def runner(idx): - kinds = ['none', 'medium', 'high'] - coup = [ 0., 0.0001, 0.001] - kind = kinds[idx] - coupling( coup[idx], kind) - data = np.genfromtxt(f'output/circ_{kind}.csv', delimiter = ',') - period_cc = np.genfromtxt(f'output/circ_dist_{kind}.csv', delimiter = ',') - data = data[:,:-1] - data = np.roll(data, -100, axis = 1) - - cells = np.zeros((2, 800,280)) - for i in tqdm.tqdm(range(0,800,1)): - Sav, rSav = run(data[i,:], period_cc[i]) - cells[:, i,:] = rSav[[0,-1],:] - np.save(f'final_cells/cells_{kind}_G1', cells) - - -# runner(0) -if __name__ == '__main__': - index = np.arange(0,3,1) - pool = Pool(3) - res = pool.map(runner, index) - pool.close() - pool.join() - - - - - -# data = np.genfromtxt('nica_output/circ_none.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_none.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(0,800,1)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_none_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_medium.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_medium.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_medium_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_high.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_high.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_high_G1', cells) - - - - - - diff --git a/Fig4_S4/CC_model/G2_coupled_system.py b/Fig4_S4/CC_model/G2_coupled_system.py deleted file mode 100644 index a2b1ac1..0000000 --- a/Fig4_S4/CC_model/G2_coupled_system.py +++ /dev/null @@ -1,200 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('K_PAPER') -from tools import * -from scipy.signal import find_peaks -import tqdm -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from poincare_oscillator import coupling -from multiprocessing import Pool -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' -''' -author @malthenielse -malthe.nielsen@nbi.ku.dk -''' - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def log2(x, xhalf, beta): - return 1/(1 + np.exp(beta*(x - xhalf))) - -def run(circ, period): - circ = np.roll(circ, np.random.randint(30,200,1)[0]) - RT = 0 - circ = (1/2 + 1/2*circ)**1.6 - # circ /= np.max(circ) - circ = np.nan_to_num(circ) - - Vol = 2 - - G1 = 3 - G2 = 0 - S = 0 - Div = 0 - AM = 0 - # RT = 0 - Sav = np.zeros((100000, 6)) - - k1 = 50 - k2 = 20 - k3 = 200 - k4 = 9 - k5 = 0.7 - k6 = 0.5 - k6 = 0.9 - # k6 = 0.01 - k7 = 1 - k8 = 5e-3 - - alpha = .3 - # alpha = 0 - - GMAX = period*3.8 + stats.expon.rvs(0, 1/.008,1)[0] - GMAX = int(GMAX) - - - - rSav = np.zeros((4,280)) - - for i in range(100000): - - ptime = int((RT/.1)) - - dG1p = k1*(1-S)*Vol - dG1n = k2 * G1 / (G1 + k3)*Vol - dS = k5*max(G1-GMAX*Vol,0)*(1-S)#*circ[ptime] - - dG2p = k4*(circ[ptime]*alpha + (1-alpha))*S*Vol - dDiv = k6 * log2(G2, 150*Vol, -1)*Vol#*circ[ptime] - dAM = (1-S)*k7*Vol - dA = 0*k8*log2(AM, 10*Vol, -1)*Vol - - rates = np.array([dG1p, dG1n, dG2p, dDiv, dS, dAM, dA]) - trates = np.sum(rates) - DT = 1/trates * np.log(1/np.random.uniform(0,1,1)[0]) - rates /= trates - - ran = np.random.uniform(0,1,1)[0] - tsum = 0 - idx = 0 - for j in range(7): - tsum += rates[j] - if ran < tsum: - idx = j - break - - if idx == 0: - G1 += 1 - - elif idx == 1: - G1 -= 1 - - elif idx == 2: - G2 += 1 - - elif idx == 3: - Div += 1 - G2 = 0 - S = 0 - rSav[0,int(RT*2):] = Div - - elif idx == 4: - S = 1 - AM = 1 - - elif idx == 5: - AM += 1 - - elif idx == 6: - Sav[i:, :] = np.nan - rSav[-1,int(RT*2)] = 1 - - break - - Sav[i,0] = G1 - Sav[i,1] = G2 - Sav[i,2] = circ[ptime] - Sav[i,3] = RT - Sav[i,4] = dA - Sav[i,5] = AM - rSav[1,int(RT*2)] = G1 - rSav[2,int(RT*2)] = G2 - - - RT += DT - - if RT > 140: - Sav[i:, :] = np.nan - break - - return Sav, rSav - -def runner(idx): - kinds = ['none', 'medium', 'high'] - kind = kinds[idx] - data = np.genfromtxt(f'output/circ_{kind}.csv', delimiter = ',') - period_cc = np.genfromtxt(f'output/circ_dist_{kind}.csv', delimiter = ',') - data = data[:,:-1] - data = np.roll(data, -100, axis = 1) - - cells = np.zeros((2, 800,280)) - for i in tqdm.tqdm(range(0,800,1)): - Sav, rSav = run(data[i,:], period_cc[i]) - cells[:, i,:] = rSav[[0,-1],:] - np.save(f'final_cells/cells_{kind}_G2', cells) - - -if __name__ == '__main__': - index = np.arange(0,3,1) - pool = Pool(3) - res = pool.map(runner, index) - pool.close() - pool.join() - - - - - -# # coupling(0.000001, 'none') -# data = np.genfromtxt('nica_output/circ_none.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_none.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(0,800,1)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_none_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_medium.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_medium.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_medium_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_high.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_high.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_high_G1', cells) - - - - - - diff --git a/Fig4_S4/CC_model/M_coupled_system.py b/Fig4_S4/CC_model/M_coupled_system.py deleted file mode 100644 index b24976a..0000000 --- a/Fig4_S4/CC_model/M_coupled_system.py +++ /dev/null @@ -1,200 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('K_PAPER') -from tools import * -from scipy.signal import find_peaks -import tqdm -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from poincare_oscillator import coupling -from multiprocessing import Pool -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' -''' -author @malthenielse -malthe.nielsen@nbi.ku.dk -''' - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def log2(x, xhalf, beta): - return 1/(1 + np.exp(beta*(x - xhalf))) - -def run(circ, period): - circ = np.roll(circ, np.random.randint(30,200,1)[0]) - RT = 0 - circ = (1/2 + 1/2*circ)**1.6 - # circ /= np.max(circ) - circ = np.nan_to_num(circ) - - Vol = 2 - - G1 = 3 - G2 = 0 - S = 0 - Div = 0 - AM = 0 - # RT = 0 - Sav = np.zeros((100000, 6)) - - k1 = 50 - k2 = 20 - k3 = 200 - k4 = 9 - k5 = 0.7 - k6 = 0.5 - k6 = 0.9 - # k6 = 0.01 - k7 = 1 - k8 = 5e-3 - - alpha = .9 - # alpha = 0 - - GMAX = period*3.8 + stats.expon.rvs(0, 1/.008,1)[0] - GMAX = int(GMAX) - - - - rSav = np.zeros((4,280)) - - for i in range(100000): - - ptime = int((RT/.1)) - - dG1p = k1*(1-S)*Vol - dG1n = k2 * G1 / (G1 + k3)*Vol - dS = k5*max(G1-GMAX*Vol,0)*(1-S) - - dG2p = k4*S*Vol - dDiv = k6 * log2(G2, 150*Vol, -1)*Vol*circ[ptime] - dAM = (1-S)*k7*Vol - dA = 0*k8*log2(AM, 10*Vol, -1)*Vol - - rates = np.array([dG1p, dG1n, dG2p, dDiv, dS, dAM, dA]) - trates = np.sum(rates) - DT = 1/trates * np.log(1/np.random.uniform(0,1,1)[0]) - rates /= trates - - ran = np.random.uniform(0,1,1)[0] - tsum = 0 - idx = 0 - for j in range(7): - tsum += rates[j] - if ran < tsum: - idx = j - break - - if idx == 0: - G1 += 1 - - elif idx == 1: - G1 -= 1 - - elif idx == 2: - G2 += 1 - - elif idx == 3: - Div += 1 - G2 = 0 - S = 0 - rSav[0,int(RT*2):] = Div - - elif idx == 4: - S = 1 - AM = 1 - - elif idx == 5: - AM += 1 - - elif idx == 6: - Sav[i:, :] = np.nan - rSav[-1,int(RT*2)] = 1 - - break - - Sav[i,0] = G1 - Sav[i,1] = G2 - Sav[i,2] = circ[ptime] - Sav[i,3] = RT - Sav[i,4] = dA - Sav[i,5] = AM - rSav[1,int(RT*2)] = G1 - rSav[2,int(RT*2)] = G2 - - - RT += DT - - if RT > 140: - Sav[i:, :] = np.nan - break - - return Sav, rSav - -def runner(idx): - kinds = ['none', 'medium', 'high'] - kind = kinds[idx] - data = np.genfromtxt(f'output/circ_{kind}.csv', delimiter = ',') - period_cc = np.genfromtxt(f'output/circ_dist_{kind}.csv', delimiter = ',') - data = data[:,:-1] - data = np.roll(data, -100, axis = 1) - - cells = np.zeros((2, 800,280)) - for i in tqdm.tqdm(range(0,800,1)): - Sav, rSav = run(data[i,:], period_cc[i]) - cells[:, i,:] = rSav[[0,-1],:] - np.save(f'final_cells/cells_{kind}_M', cells) - - -if __name__ == '__main__': - index = np.arange(0,3,1) - pool = Pool(3) - res = pool.map(runner, index) - pool.close() - pool.join() - - - - - -# # coupling(0.000001, 'none') -# data = np.genfromtxt('nica_output/circ_none.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_none.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(0,800,1)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_none_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_medium.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_medium.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_medium_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_high.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_high.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_high_G1', cells) - - - - - - diff --git a/Fig4_S4/CC_model/S_coupled_system.py b/Fig4_S4/CC_model/S_coupled_system.py deleted file mode 100644 index bfdd8bd..0000000 --- a/Fig4_S4/CC_model/S_coupled_system.py +++ /dev/null @@ -1,201 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('K_PAPER') -from tools import * -from scipy.signal import find_peaks -import tqdm -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from poincare_oscillator import coupling -from multiprocessing import Pool -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' -''' -author @malthenielse -malthe.nielsen@nbi.ku.dk -''' - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def log2(x, xhalf, beta): - return 1/(1 + np.exp(beta*(x - xhalf))) - -def run(circ, period): - circ = np.roll(circ, np.random.randint(30,200,1)[0]) - RT = 0 - circ = (1/2 + 1/2*circ)**1.6 - # circ /= np.max(circ) - circ = np.nan_to_num(circ) - - Vol = 2 - - G1 = 3 - G2 = 0 - S = 0 - Div = 0 - AM = 0 - # RT = 0 - Sav = np.zeros((100000, 6)) - - k1 = 50 - k2 = 20 - k3 = 200 - k4 = 9 - k5 = 0.7 - k6 = 0.5 - k6 = 0.9 - # k6 = 0.01 - k7 = 1 - k8 = 5e-3 - - alpha = 0 - # alpha = 0 - - GMAX = period*3.8 + stats.expon.rvs(0, 1/.005,1)[0] - GMAX = int(GMAX) - - - - rSav = np.zeros((4,280)) - - for i in range(100000): - - ptime = int((RT/.1)) - - dG1p = k1*(1-S)*Vol - dG1n = k2 * G1 / (G1 + k3)*Vol - dS = k5*max(G1-GMAX*Vol,0)*(1-S)*circ[ptime] - - dG2p = k4*S*Vol - dDiv = k6 * log2(G2, 150*Vol, -1)*Vol#*circ[ptime] - dAM = (1-S)*k7*Vol - dA = 0*k8*log2(AM, 10*Vol, -1)*Vol - - rates = np.array([dG1p, dG1n, dG2p, dDiv, dS, dAM, dA]) - trates = np.sum(rates) - DT = 1/trates * np.log(1/np.random.uniform(0,1,1)[0]) - rates /= trates - - ran = np.random.uniform(0,1,1)[0] - tsum = 0 - idx = 0 - for j in range(7): - tsum += rates[j] - if ran < tsum: - idx = j - break - - if idx == 0: - G1 += 1 - - elif idx == 1: - G1 -= 1 - - elif idx == 2: - G2 += 1 - - elif idx == 3: - Div += 1 - G2 = 0 - S = 0 - rSav[0,int(RT*2):] = Div - - elif idx == 4: - S = 1 - AM = 1 - - elif idx == 5: - AM += 1 - - elif idx == 6: - Sav[i:, :] = np.nan - rSav[-1,int(RT*2)] = 1 - - break - - Sav[i,0] = G1 - Sav[i,1] = G2 - Sav[i,2] = circ[ptime] - Sav[i,3] = RT - Sav[i,4] = dA - Sav[i,5] = AM - rSav[1,int(RT*2)] = G1 - rSav[2,int(RT*2)] = G2 - - - RT += DT - - if RT > 140: - Sav[i:, :] = np.nan - break - - return Sav, rSav - -def runner(idx): - kinds = ['none', 'medium', 'high'] - kind = kinds[idx] - print(kind) - data = np.genfromtxt(f'output/circ_{kind}.csv', delimiter = ',') - period_cc = np.genfromtxt(f'output/circ_dist_{kind}.csv', delimiter = ',') - data = data[:,:-1] - # data = np.roll(data, -100, axis = 1) - - cells = np.zeros((2, 800,280)) - for i in tqdm.tqdm(range(0,800,1)): - Sav, rSav = run(data[i,:], period_cc[i]) - cells[:, i,:] = rSav[[0,-1],:] - np.save(f'final_cells/cells_{kind}_S', cells) - - -if __name__ == '__main__': - index = np.arange(0,3,1) - pool = Pool(3) - res = pool.map(runner, index) - pool.close() - pool.join() - - - - - -# # coupling(0.000001, 'none') -# data = np.genfromtxt('nica_output/circ_none.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_none.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(0,800,1)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_none_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_medium.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_medium.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_medium_G1', cells) -# -# # coupling(0.005, 'high_alt') -# data = np.genfromtxt('nica_output/circ_high.csv', delimiter = ',') -# period_cc = np.genfromtxt('nica_output/circ_dist_high.csv', delimiter = ',') -# data = data[:,:-1] -# data = np.roll(data, -100, axis = 1) -# cells = np.zeros((2, 800,280)) -# for i in tqdm.tqdm(range(800)): -# Sav, rSav = run(data[i,:], period_cc[i]) -# cells[:, i,:] = rSav[[0,-1],:] -# np.save('final_cells/cells_high_G1', cells) - - - - - - diff --git a/Fig4_S4/CC_model/__pycache__/poincare_oscillator.cpython-311.pyc b/Fig4_S4/CC_model/__pycache__/poincare_oscillator.cpython-311.pyc deleted file mode 100644 index aadea3d..0000000 Binary files a/Fig4_S4/CC_model/__pycache__/poincare_oscillator.cpython-311.pyc and /dev/null differ diff --git a/Fig4_S4/CC_model/__pycache__/tools.cpython-311.pyc b/Fig4_S4/CC_model/__pycache__/tools.cpython-311.pyc deleted file mode 100644 index 18028b0..0000000 Binary files a/Fig4_S4/CC_model/__pycache__/tools.cpython-311.pyc and /dev/null differ diff --git a/Fig4_S4/CC_model/cells_plotter.py b/Fig4_S4/CC_model/cells_plotter.py deleted file mode 100644 index 6f55a65..0000000 --- a/Fig4_S4/CC_model/cells_plotter.py +++ /dev/null @@ -1,277 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -plt.style.use('K_PAPER') -from tools import * -from appstatpy.ExternalFunctions import * -from scipy import stats, special -from iminuit import Minuit -from scipy.optimize import curve_fit -# output = '/home/nordentoft/Nextcloud/Manuscripts/Circadian-CellCycle/Malthe_figures/' -from matplotlib.colors import LinearSegmentedColormap - -def create_custom_colormap(num_colors): #input: maximum number divisions - all_colors = [ - (135/255, 206/255, 250/255), # skyblue - (255/255, 228/255, 181/255), # moccasin - (216/255, 191/255, 216/255), # thistle - (240/255, 128/255, 128/255), # lightcoral - (95/255, 158/255, 160/255), # cadetblue - (255/255, 215/255, 0/255), # gold - ] - - if num_colors not in range(3, 7): - raise ValueError("Number of colors must be between 3 and 6 inclusive") - - colors = all_colors[:num_colors] - - cmap_name = 'custom_colormap' - return LinearSegmentedColormap.from_list(cmap_name, colors, N=256) - -ccmap = create_custom_colormap(6) - -def exponnorm(x, mu, sigma, lam, N): - erf = 1 - special.erf((mu + lam*sigma**2 - x)/(np.sqrt(2)*sigma)) - return N * lam/2*np.exp(lam/2*(2*mu + lam*sigma**2 - 2*x))*erf - -def lin_func(x, a, b): - return a*x + b - -# kind = 'S' -kind = 'G1' -# kind = 'G2' -# kind = 'M' - -ch = np.load(f'final_cells/cells_high_{kind}.npy') -cn = np.load(f'final_cells/cells_none_{kind}.npy') -cm = np.load(f'final_cells/cells_medium_{kind}.npy') - -dch = ch[1,:,:] -ch = ch[0,:,:] -dcn = cn[1,:,:] -cn = cn[0,:,:] -dcm = cm[1,:,:] -cm = cm[0,:,:] - -def div_matrix(mat, string): - dmat = np.diff(mat, axis = 1) - dmat = dmat[:,40:] - np.savetxt('dmats/'+string+'.csv', dmat, delimiter = ',') -# -div_matrix(ch, 'high' + '_' + kind ) -div_matrix(cn, 'none' + '_' + kind) -div_matrix(cm, 'midle' + '_' + kind) - - -# exit() - -# ch = np.load('./nica_output/pmat_high.npy') -# cn = np.load('./nica_output/pmat_none.npy') -# ch = np.cumsum(ch, axis = 1) -# cn = np.cumsum(cn, axis = 1) - - -# dch = np.cumsum(np.sum(dch, axis =0)) -# dcn = np.cumsum(np.sum(dcn, axis =0)) - - -# time = np.arange(20.5, ch.shape[1]//2 - 40, .5) - - -# def sin_lin(time, alpha, beta, A0, T, phase): -# return alpha*time + beta + A0*np.sin(2*np.pi/T*time + phase) -# -# pm = prolif_matrix(ch)/400 -# pm = pm[55:] -# time = np.arange(0, len(pm)/2, .5) -# fig, ax = plt.subplots(1,2, figsize = (14,7)) -# ax[0].plot(time, pm, color = 'black') -# p0, _ = curve_fit(sin_lin, time, pm, p0 = [1/24, 0, 10, 20, 0]) -# ax[0].plot(time, sin_lin(time, *p0), color = 'black', ls = '--') -# print(p0, 1/p0[0]) -# -# ax[1].plot(time, pm - (p0[0]*time + p0[1])) -# ax[1].plot(time, p0[2]*np.sin(2*np.pi/p0[3]*time + p0[4])) -# -# plt.show() -# exit() - -# fig, ax = plt.subplots(1,1, figsize = (7,7)) - - - - -idx = master_sort(ch) -ch = ch[idx, :] - -idx = master_sort(cn) -cn = cn[idx, :] - -idx = master_sort(cm) -cm = cm[idx, :] - -fig, ax = plt.subplots(1,3, figsize = (7,10), sharex = True, sharey = True) -ax[0].imshow(ch, aspect = 'auto', extent = [-20,120, 0, 400], cmap = ccmap, vmin = 0, vmax = 5) -ax[1].imshow(cm, aspect = 'auto',extent = [-20,120, 0, 400] , cmap = ccmap, vmin = 0, vmax = 5) -ax[2].imshow(cn, aspect = 'auto',extent = [-20,120, 0, 400] , cmap = ccmap, vmin = 0, vmax = 5) -ax[0].set_xlim(0,120) -ax[0].set(title = 'Untreated', xlabel = 'Time [hr]', ylabel = 'Cell #') -ax[1].set(title = '5uM', xlabel = 'Time [hr]') -ax[2].set(title = '10uM', xlabel = 'Time [hr]') -fig.savefig(output + 'prolif_matrix_simulation.svg', dpi = 500) - -# plt.show() - -time = np.arange(.5, ch.shape[1]/2, .5) -# print(ch.shape) - - -fig, ax = plt.subplots(1,1, figsize = (7,5)) -pch = prolif_matrix(ch) -pcm = prolif_matrix(cm) -pcn = prolif_matrix(cn) -max_val = max([pch[-1], pcm[-1], pcn[-1]]) -max_val = ch.shape[0] - -grate = [] - -ax.plot(time[40:], pch[40:]/max_val, color = 'black', label = 'High') -p0, _ = curve_fit(lin_func, time[40:],prolif_matrix(ch)[40:]/max_val) -p00, _ = curve_fit(lin_func, time[40:],prolif_matrix(ch)[40:]/800) -grate.append(p0[0]) -# print(time[40]) -# print(p0) -# ax.plot(time[40:], lin_func(time[40:], *p0), color = 'black', ls = '--', label = f'Slope = {np.round(1/p00[0],1)}') - -ax.plot(time[40:], pcm[40:]/max_val, color = 'red', label = 'low') -p0, _ = curve_fit(lin_func, time[40:], prolif_matrix(cm)[40:]/max_val) -p00, _ = curve_fit(lin_func, time[40:], prolif_matrix(cm)[40:]/800) -grate.append(p0[0]) -# ax.plot(time[40:], lin_func(time[40:], *p0), color = 'red', ls = '--', label = f'Slope = {np.round(1/p00[0],1)}') - -ax.plot(time[40:], pcn[40:]/max_val, color = 'green', label = 'None') -p0, _ = curve_fit(lin_func, time[130:],prolif_matrix(cn)[130:]/max_val) -p00, _ = curve_fit(lin_func, time[130:],prolif_matrix(cn)[130:]/800) -grate.append(p0[0]) -# print(1/p0) -# ax.plot(time[40:], lin_func(time[40:], *p0), color = 'green', ls = '--', label = f'Slope = {np.round(1/p00[0],1)}') - -ax.set_xticks(np.linspace(20,150,6), [0, 24, 48, 72, 96, 120]) - -ax.legend() - -# fig, ax = plt.subplots(1,1, figsize = (7,7)) -# ax[1].plot(time[14:-5], run_mean(np.diff(run_mean(prolif_matrix(ch),10)), 10), color = 'black') -# ax[1].plot(time[14:-5], run_mean(np.diff(run_mean(prolif_matrix(cn),10)), 10), color = 'green') -# -# ax[1].set(xlabel = 'Time [hr]', ylabel = 'Signal[AU]', title = 'Detrended growth') -ax.set(xlabel = 'Time [hr]', ylabel = 'Total fate', title = 'Growth') -fig.savefig(output + 'Growth_simulation.svg', dpi = 500) - - - -fit = {'ExpNorm' : ['$\mu$', '$\sigma_{\mu}$']} -fig, ax = plt.subplots(1,1, figsize = (7,7)) -bins = np.arange(10, 100,1) - -err_bin = [] - -res = np.zeros((3,3)) - -lib, _ = extract_cell_division(ch) -# lib['prolif_total'] = np.delete(np.where(lib['prolif_total'] > 70)[0],lib['prolif_total']) -lib['prolif_total'] = lib['prolif_total'][np.where(lib['prolif_total'] < 90)[0]] - -# print(lib['prolif_total']) -# ax.hist(lib['prolif_total']/2, bins = bins, histtype = 'step', lw = 2, density = True, color = 'black') -ax.hist(np.array(lib['prolif_total'])/2, histtype = 'step', color = 'black', bins = np.arange(10,60,.5), density = True) -freq, bins = np.histogram(lib['prolif_total'], bins = np.arange(10,60,.5), density = False) -freq_err = np.sqrt(freq) -bin_center = (bins[:-1] + bins[1:])/4 -bw = bin_center[1] - bin_center[0] -chi_im = Chi2Regression(exponnorm, bin_center[freq > 0], freq[freq > 0], freq_err[freq > 0]) -chi_im = Minuit(chi_im, mu = 20, sigma = 2.1, lam = 0.1, N = np.sum(freq)*bw) -chi_im.errordef = 1 -chi_im.migrad() -ullh_x = np.linspace(0,100,1000) -print(chi_im.values[:]) -ax.plot(ullh_x, exponnorm(ullh_x, *chi_im.values[:-1], N = 1), color = 'black') -mu_true = chi_im.values[0] + 1/chi_im.values[2] -chi2 = stats.chi2.sf(chi_im.fval, len(freq[freq >0])) -err = np.sqrt(chi_im.errors[0]**2 + (1/chi_im.values[2]**2)*chi_im.errors[2]**2) -print(mu_true, err, chi2) -res[0,:] = chi_im.values[:3] -err_bin.append(err) - -fit['Untreated'] = [mu_true, err] - -lib, _ = extract_cell_division(cm) -lib['prolif_total'] = lib['prolif_total'][np.where(lib['prolif_total'] < 90)[0]] -ax.hist(np.array(lib['prolif_total'])/2, histtype = 'step', color = 'red', bins = np.arange(10,60,.5), density = True) -freq, bins = np.histogram(lib['prolif_total'], bins = np.arange(10,60,.5), density = False) -freq_err = np.sqrt(freq) -bin_center = (bins[:-1] + bins[1:])/4 -bw = bin_center[1] - bin_center[0] -chi_im = Chi2Regression(exponnorm, bin_center[freq > 0], freq[freq > 0], freq_err[freq > 0]) -chi_im = Minuit(chi_im, mu = 20, sigma = 2.1, lam = 0.1, N = np.sum(freq)*bw) -chi_im.errordef = 1 -chi_im.migrad() -ullh_x = np.linspace(0,100,1000) -print(chi_im.values[:]) -ax.plot(ullh_x, exponnorm(ullh_x, *chi_im.values[:-1], N = 1), color = 'red') -mu_true = chi_im.values[0] + 1/chi_im.values[2] -chi2 = stats.chi2.sf(chi_im.fval, len(freq[freq >0])) -err = np.sqrt(chi_im.errors[0]**2 + (1/chi_im.values[2]**2)*chi_im.errors[2]**2) -print(mu_true, err, chi2) -res[1,:] = chi_im.values[:3] -fit['5uM'] = [mu_true, err] -err_bin.append(err) - -lib, _ = extract_cell_division(cn) -lib['prolif_total'] = lib['prolif_total'][np.where(lib['prolif_total'] < 90)[0]] -ax.hist(np.array(lib['prolif_total'])/2, histtype = 'step', color = 'green', bins = np.arange(10,60,.5), density = True) -freq, bins = np.histogram(lib['prolif_total'], bins = np.arange(10,60,.5), density = False) -freq_err = np.sqrt(freq) -bin_center = (bins[:-1] + bins[1:])/4 -bw = bin_center[1] - bin_center[0] -chi_im = Chi2Regression(exponnorm, bin_center[freq > 0], freq[freq > 0], freq_err[freq > 0]) -chi_im = Minuit(chi_im, mu = 18, sigma = 2.3, lam = 0.1, N = np.sum(freq)*bw) -chi_im.errordef = 1 -chi_im.migrad() -ullh_x = np.linspace(0,100,1000) -print(chi_im.values[:]) -ax.plot(ullh_x, exponnorm(ullh_x, *chi_im.values[:-1], N = 1), color = 'green') -mu_true = chi_im.values[0] + 1/chi_im.values[2] -chi2 = stats.chi2.sf(chi_im.fval, len(freq[freq >0])) -err = np.sqrt(chi_im.errors[0]**2 + (1/chi_im.values[2]**2)*chi_im.errors[2]**2) -print(mu_true, err, chi2) -err_bin.append(err) - -res[2,:] = chi_im.values[:3] -fit['10uM'] = [mu_true, err] -# ax.plot(ullh_x, exponnorm(ullh_x, 24, 5, .1, N = 1), color = 'red') - -grate = np.array([grate]) - -res = np.concatenate((res, grate), axis = 0) -np.save(f'res_{kind}.npy', res) - -text = nice_string_output(fit, extra_spacing=0, decimals=3) -add_text_to_ax(.5, .5, text, ax, fontsize = 11) -ax.set(xlabel = 'Period [hr]', ylabel = 'p', title = 'Intermiotic time', xlim = (15, 55)) -# fig.savefig(output + 'IM_historgram_simulation.svg', dpi = 500) - -# time = np.arange(0, len(avg)/2, .5) -# axx[0].plot(time, avg/max_val, color = colors[i], alpha = 1) -# axx[0].plot(time, avg/351, color = colors[i], alpha = 1) - -# p0, _ = curve_fit(lin_func, time[:120], avg[:120]/max_val, p0 = [1, 1]) -# p0, _ = curve_fit(lin_func, time[:90], avg[:90]/351, p0 = [1, 1]) - - -plt.show() - - -print(err_bin) - - - diff --git a/Fig4_S4/CC_model/poincare_oscillator.py b/Fig4_S4/CC_model/poincare_oscillator.py deleted file mode 100644 index d0b2418..0000000 --- a/Fig4_S4/CC_model/poincare_oscillator.py +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 14 10:16:45 2023 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -from scipy import stats -import tqdm -# import seaborn as sns -# from scipy.signal import find_peaks - -np.random.seed(42) - - -def coupling(kappa, kind): -#Parameters - Nosc = 800 - dt = 0.1 - tf = 150 - - gg = 1 - ll = 1 - a0 = 1 - acc = 1 - -#Coupling parameters - # kappa = 0.0001 #extracellular circadian coupling - eps = 0.005 #intracellular circadian - cell cycle - -#Vectors initialization - t = np.arange(0,tf,step=dt) - X = np.zeros((Nosc,len(t) + 1)) - XX = np.zeros((Nosc,len(t))) - Y = np.zeros((Nosc,len(t))) - YY = np.zeros((Nosc,len(t))) - -#Initial position - # X[:,0] = np.random.uniform(.0,.3,size=(Nosc)) - # XX[:,0] = np.random.uniform(-1,1,size=(Nosc)) - # Y[:,0] = np.random.uniform(-.0,.3,size=(Nosc)) - # YY[:,0] = np.random.uniform(-1,1,size=(Nosc)) - X[:,0] = np.random.uniform(.0,.1,size=(Nosc)) - XX[:,0] = np.random.uniform(-1,1,size=(Nosc)) - Y[:,0] = np.random.uniform(-.0,.1,size=(Nosc)) - YY[:,0] = np.random.uniform(-1,1,size=(Nosc)) - -#Period distribution - # mu1, sigma1 = 23, np.sqrt(6) - mu1, sigma1 = 22, 3 - period_circ = np.random.normal(mu1, sigma1, Nosc)#+ stats.expon.rvs(0, 1/.1,Nosc) - mu2, sigma2 = 26, np.sqrt(6) - # mu2, sigma2 = 26, 2.5 - period_cell = np.random.normal(mu2, sigma2, Nosc)#+ stats.expon.rvs(0, 1/.1,Nosc) - -# combined_data = np.concatenate((period_circ, period_cell)) -# bin_edges = np.histogram_bin_edges(combined_data, bins=25) - -# plt.figure(figsize=(12,10)) -# num_bins = 25 -# sns.histplot(period_circ, stat='density', kde=True, bins=bin_edges, color='darkblue', label='Circadian') -# sns.histplot(period_cell, stat='density', kde=True, bins=bin_edges, color='darkgreen', label='Cell cycle') -# plt.xlabel('Period') -# plt.legend(loc='best') -# plt.show() - -# Integration over time and population - SUMS = np.zeros((2,1500)) - for i in tqdm.tqdm(range(len(t) - 1)): - sum_x = 0 - sum_y = 0 - - for m in range(Nosc): - middle = -gg*(np.sqrt(X[m,i]**2+Y[m,i]**2)-a0) - X[m,i+1] = X[m,i]+dt*(middle*X[m,i]-2*np.pi*Y[m,i]/period_circ[m]+sum_x*kappa/(2*Nosc)) - Y[m,i+1] = Y[m,i]+dt*(middle*Y[m,i]+2*np.pi*X[m,i]/period_circ[m]+sum_y*kappa/(2*Nosc)) - - mid = -ll*(np.sqrt(XX[m,i]**2+YY[m,i]**2)-acc) - XX[m,i+1] = XX[m,i]+dt*(mid*XX[m,i]-2*np.pi*YY[m,i]/period_cell[m]+eps*(X[m,i]+XX[m,i])/2) - YY[m,i+1] = YY[m,i]+dt*(mid*YY[m,i]+2*np.pi*XX[m,i]/period_cell[m]+eps*(Y[m,i]+YY[m,i])/2) - - for k in range(Nosc): - sum_x = sum_x+X[k,i] - sum_y = sum_y+Y[k,i] - - SUMS[0,i] = sum_x - SUMS[1,i] = sum_y - - - fig, ax = plt.subplots(1,1, figsize = (8,6)) - for i in range(800): - X[i,:] = X[i,:]/max(X[i,:]) - # X[i,:] = np.roll(X[i,:], - np.random.randint(30, 100,1)) - # inds, _ = find_peaks(X[i,:]) - # print(np.diff(inds)*dt, period_cell[i], period_circ[i]) - # X[i,-1] = period_circ[i] - # ax.plot(X[i,:], color = 'black', alpha = .08) - - - -# ax.plot(SUMS[0,:]/np.max(abs(SUMS[0,:])), color = 'red') -# np.savetxt('half_coupled.csv', X/np.max(X), delimiter = ',') -# np.savetxt('coupled.csv', X, delimiter = ',') -# plt.show() - np.savetxt(f'output/circ_{kind}.csv', X, delimiter = ',') - np.savetxt(f'output/circ_dist_{kind}.csv', period_circ, delimiter = ',') - diff --git a/Fig4_S4/CC_model/tools.py b/Fig4_S4/CC_model/tools.py deleted file mode 100644 index f1583d5..0000000 --- a/Fig4_S4/CC_model/tools.py +++ /dev/null @@ -1,106 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt -# plt.style.use('science') -from scipy import stats - -def master_sort(CL): - N = CL.shape[0] - CL_max = CL.max(axis = 1) - pro_set = np.array(list(set(CL_max))).astype(int) - key = [] - for pro in pro_set: - idx = np.where(CL_max == pro)[0] - if pro == 0: - # FP_index = idx.copy() - # key = list(idx) - key.append(idx) - continue - - FP_index = [] - for index in idx: - FP = np.nonzero(CL[index,:])[0][0] - FP_index.append(FP) - FP_index = np.array(FP_index) - FP_key = np.argsort(FP_index) - FP_index = FP_index[FP_key] - idx = idx[FP_key] - key.append(idx) - key = np.hstack(key) - return key - -def extract_cell_division(CL): - N = CL.shape[0] - CL_max = CL.max(axis = 1) - pro_set = np.array(list(set(CL_max))).astype(int) - lib = {} - lib_time = {} - total_peaks = [] - total_times = [] - for pro in pro_set: - idx = np.where(CL_max == pro)[0] - if pro < 2: - continue - peak_list = [] - time_list = [] - for index in idx: - dCL = np.diff(CL[index,:]) - peaks = np.where(dCL == 1)[0] - dPeaks = np.diff(peaks) - peak_list.append(dPeaks) - time_list.append(peaks[:-1]) - peak_list = np.hstack(peak_list) - time_list = np.hstack(time_list) - total_peaks.append(peak_list) - total_times.append(time_list) - lib[f"prolif_{pro}"] = peak_list - lib_time[f"prolif_{pro}"] = time_list - total_peaks = np.hstack(total_peaks) - total_times = np.hstack(total_times) - lib["prolif_total"] = total_peaks - lib_time["prolif_total"] = total_times - - return lib, lib_time - -def cells_correlation(CL): - period_lib = {} - times_lib = {} - CL_max = CL.max(axis = 1) - # print(CL_max) - for i in range(1,int(np.max(CL_max)),1): - period_lib[f'div_{i}'] = [] - times_lib[f'div_{i}'] = [] - - # print(CL.shape) - - for i in range(CL.shape[0]): - if CL_max[i] < 2: - continue - diff = np.diff(CL[i,:]) - idx = np.where(diff == 1)[0] - for j in range(len(idx)-1): - period = idx[j+1] - idx[j] - period_lib[f'div_{j+1}'].append(period) - times_lib[f'div_{j+1}'].append(idx[j]) - # for i in range(1,int(np.max(CL_max)),1): - # if len(times_lib[f'div_{i}']) < 1: - # continue - # print(stats.pearsonr(times_lib[f'div_{i}'], period_lib[f'div_{i}'])[0]) - return times_lib, period_lib - - - - - - - - -def prolif_matrix(CL): - dCL = np.diff(CL, axis = 1) - sumCL = np.sum(dCL, axis = 0) - return np.cumsum(sumCL) - -def run_mean(x, N): - return np.convolve(x, np.ones(N)/N, mode='valid') - - - diff --git a/Fig4_S4/Cum_distr_osc_periods.py b/Fig4_S4/Cum_distr_osc_periods.py deleted file mode 100644 index 1dd6073..0000000 --- a/Fig4_S4/Cum_distr_osc_periods.py +++ /dev/null @@ -1,283 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:59:26 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -import random -from pyboat import WAnalyzer -import seaborn as sns -from scipy.optimize import curve_fit - -def exponential_func(t, N0, r): - K = max(df[j]) - return K/(1+((K-N0)/N0)*np.exp(-r*t)) - -def linear_func(t, a, b): - return a*t+b - -dt = 0.5 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -path = '.../Division_matrix/' -output = '.../Fig4/' -dose = ['untreated','5uM','10uM'] -density = ['high'] -channel2 = 'cell_cycle' - -#%%One density -for i in density: - - df = pd.DataFrame() - fit_params = [] - cov_error = [] - - fit_params_l = [] - cov_error_l = [] - - plt.figure(figsize=(12,10)) - - for j in dose: - condition = str(i)+'_density_'+str(j)+'_' - print(condition) - - #Division matrix - data = pd.read_excel(path+'divisions_'+str(condition)+r'.xlsx') - time = data.index.values - print(len(data.columns)) - - median_distr = pd.DataFrame() - for a in range(100): - num_Ids = random.sample(list(data.columns), 750) #high=750, medium=112, low=49 - time_points = len(time) - division_profile = np.zeros((len(num_Ids),time_points)) - - division_profile = np.array([data[col].values for col in num_Ids]) - - distr = (np.sum(division_profile, axis=0)) - median_distr[a] = np.cumsum(distr) - - average_distr = median_distr.mean(axis=1) - - df[j] = average_distr - maximum = (df.values.max()) - - params_l, covariance_l = curve_fit(linear_func, time[0:150], df[j][0:150]) - errors_l = np.sqrt(np.diag(covariance_l)) - y_fit_l = linear_func(time[0:150], *params_l) - fit_params_l.append(params_l[0]) - cov_error_l.append(errors_l[0]) - - params, covariance = curve_fit(exponential_func, time, df[j]) - errors = np.sqrt(np.diag(covariance)) - y_fit = exponential_func(time, *params) - fit_params.append(params[1]) - cov_error.append(errors[1]) - print(errors[1]) - - plt.plot(time*0.5, df[j]/maximum, label = str(j), linewidth=4) - plt.plot(time[24:len(time)]*0.5, y_fit[24:len(time)]/maximum, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5) - # plt.plot(time[0:150]*0.5, y_fit_l, '--', label='linear fit: {:.3f}'.format(params_l[0]), linewidth=5, alpha=0.5) - plt.xticks([0,24,48,72,96,120]) - plt.xlabel('Time(h)') - plt.ylabel('Cumulative distribution of division events') - plt.legend(loc = 'best') - # plt.savefig(output+'Cumul_distr_all_divisions_0_5_10uM.svg') - plt.show() - - print(fit_params) - print(cov_error) - - plt.figure(figsize=(10,8)) - plt.errorbar(dose, fit_params, yerr=cov_error, fmt='o') - plt.xlabel('Inhibitor') - plt.ylabel('Rate of the logistic function') - # plt.savefig(output+'Logistic_fit_growth_curves.svg') - plt.show() - - # plt.figure(figsize=(10,8)) - # plt.errorbar(dose, fit_params_l, yerr=cov_error_l, fmt='o') - # plt.xlabel('Inhibitor') - # plt.ylabel('Growth rate') - # # plt.savefig(output+'Growth_rates_curves.svg') - # plt.show() - - plt.figure(figsize=(12,10)) - for col in df: - plt.plot(time*0.5, wAn.sinc_smooth(wAn.sinc_detrend(df[col]/maximum, T_c=40), T_c=10), label = str(col), linewidth=4) - plt.xlabel('Time(h)') - plt.ylabel('Detrended cumulative distribution of division events') - plt.xticks([0,24,48,72,96,120]) - plt.legend(loc = 'best') - # plt.savefig(output+'Detrended_Cumul_distr_all_divisions_0_5_10uM.svg') - plt.show() - - periods = {} - - for col in df: - signal = wAn.sinc_detrend(df[col]/maximum, T_c=42) - wAn.compute_spectrum(signal, do_plot=False, draw_coi=False) - rd = wAn.get_maxRidge(power_thresh = 0, smoothing_wsize=20) - # wAn.draw_Ridge() - periods[col] = rd['periods'] - - plt.figure(figsize=(12,10)) - for col in df: - plt.plot(periods[col].index.values*0.5, periods[col]/periods['untreated'], label = str(col), linewidth=4) - plt.xlim([-5,65]) - plt.ylabel('Period detrended cumulative division events [hours]') - plt.legend(loc = 'best') - # plt.savefig(output+'Period_detrended_Cumul_distr_all_divisions_0_5_10uM.svg') - plt.show() - - periods_df = pd.DataFrame().from_dict(periods, orient='index') - periods_df = periods_df.T - - for col in periods_df: - print(col, np.mean(periods_df[col]), np.std(periods_df[col])) - - plt.figure(figsize=(12,10)) - sns.boxplot(data=periods_df, showfliers=False) - # for col in df: - # plt.plot(periods[col].index.values*0.5, periods[col], label = str(col), linewidth=4) - # plt.ylim([22,35]) - plt.ylabel('Period detrended cumulative division events [hours]') - # plt.legend(loc = 'best') - # plt.savefig(output+'Boxplot_period_detrended_cumul_distr_cut_yaxos.svg') - plt.show() - - -#%%Density - -density = ['high','medium','low'] - -plt.figure(figsize=(12,10)) - -df = pd.DataFrame() -fit_params = [] -cov_error = [] - -fit_params_l = [] -cov_error_l = [] - -for j in density: - condition = str(j)+'_density_untreated_' - print(condition) - - #Division matrix - data = pd.read_excel(path+'divisions_'+str(condition)+r'.xlsx') - time = data.index.values - print(len(data.columns)) - - median_distr = pd.DataFrame() - for a in range(100): - num_Ids = random.sample(list(data.columns), 89) #high=750, medium=112, low=49 - time_points = len(time) - division_profile = np.zeros((len(num_Ids),time_points)) - - division_profile = np.array([data[col].values for col in num_Ids]) - - distr = (np.sum(division_profile, axis=0)) - median_distr[a] = np.cumsum(distr) - - average_distr = median_distr.mean(axis=1) - - df[j] = average_distr - maximum = (df.values.max()) - - params_l, covariance_l = curve_fit(linear_func, time[0:150], df[j][0:150]) - errors_l = np.sqrt(np.diag(covariance_l)) - y_fit_l = linear_func(time[0:150], *params_l) - fit_params_l.append(params_l[0]) - cov_error_l.append(errors_l[0]) - - params, covariance = curve_fit(exponential_func, time, df[j]) - errors = np.sqrt(np.diag(covariance)) - y_fit = exponential_func(time, *params) - fit_params.append(params[1]) - cov_error.append(errors[1]) - print(errors[1]) - - - plt.plot(time*0.5, df[j]/maximum, label = str(j), linewidth=4) - plt.plot(time[24:len(time)]*0.5, y_fit[24:len(time)]/maximum, '--', label='exp fit: {:.3f}'.format(params[1]), linewidth=5, alpha=0.5) - # plt.plot(time[0:150]*0.5, y_fit_l, '--', label='linear fit: {:.3f}'.format(params_l[0]), linewidth=5, alpha=0.5) - -print(fit_params) -print(cov_error) - -plt.xticks([0,24,48,72,96,120]) -plt.xlabel('Time(h)') -plt.ylabel('Cumulative distribution of division events') -plt.legend(loc = 'best') -# plt.savefig(output+'Cumul_distr_all_divisions_change_density.svg') -plt.show() - -plt.figure(figsize=(10,8)) -plt.errorbar(dose, fit_params, yerr=cov_error, fmt='o') -plt.xlabel('Inhibitor') -plt.ylabel('Rate of the logistic function') -# plt.savefig(output+'Logistic_fit_growth_curves_density.svg') -plt.show() - -# plt.figure(figsize=(10,8)) -# plt.errorbar(dose, fit_params_l, yerr=cov_error_l, fmt='o') -# plt.xlabel('Inhibitor') -# plt.ylabel('Growth rate') -# # plt.savefig(output+'Growth_rates_curves.svg') -# plt.show() - -plt.figure(figsize=(12,10)) -for col in df: - plt.plot(time*0.5, wAn.sinc_smooth(wAn.sinc_detrend(df[col]/maximum, T_c=40), T_c=10), label = str(col), linewidth=4) -plt.xlabel('Time(h)') -plt.ylabel('Detrended cumulative distribution of division events') -plt.xticks([0,24,48,72,96,120]) -plt.legend(loc = 'best') -# plt.savefig(output+'Detrended_Cumul_distr_change_density.svg') -plt.show() - -periods = {} - -for col in df: - signal = wAn.sinc_detrend(df[col]/maximum, T_c=42) - wAn.compute_spectrum(signal, do_plot=False) - rd = wAn.get_maxRidge(power_thresh = 0, smoothing_wsize=20) - # wAn.draw_Ridge() - periods[col] = rd['periods'] - -# plt.figure(figsize=(12,10)) -# for col in df: -# plt.plot(periods[col].index.values*0.5, periods[col]/periods['untreated'], label = str(col), linewidth=4) -# plt.xlim([-5,65]) -# plt.ylabel('Period detrended cumulative division events [hours]') -# plt.legend(loc = 'best') -# # plt.savefig(output+'Period_detrended_Cumul_distr_all_divisions_0_5_10uM.svg') -# plt.show() - -periods_df = pd.DataFrame().from_dict(periods, orient='index') -periods_df = periods_df.T - -for col in periods_df: - print(col, np.mean(periods_df[col]), np.std(periods_df[col])) - -plt.figure(figsize=(12,10)) -sns.boxplot(data=periods_df, showfliers=False) -# for col in df: -# plt.plot(periods[col].index.values*0.5, periods[col], label = str(col), linewidth=4) -# plt.ylim([22,35]) -plt.ylabel('Period detrended cumulative division events [hours]') -# plt.legend(loc = 'best') -# plt.savefig(output+'Boxplot_period_detrended_cumul_distr_change_density.svg') -plt.show() \ No newline at end of file diff --git a/Fig4_S4/IMT_distributions.py b/Fig4_S4/IMT_distributions.py deleted file mode 100644 index b2e3c60..0000000 --- a/Fig4_S4/IMT_distributions.py +++ /dev/null @@ -1,140 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:57:48 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -import seaborn as sns - -dt = 0.5 # the sampling interval, 0.5hours -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -path = '.../Division_matrix/' -path2 = '.../Raw_data/' -output = '.../Fig4/' -dose = ['untreated', '5uM', '10uM'] -density = ['high'] -channel1 = 'circadian' - -colors = {'untreated':'gold', '5uM':'tab:orange', '10uM':'tab:blue'} - -#%% Changing inhibitor -imt_all_dose = {} - -for i in density: - df = pd.DataFrame() - - for j in dose: - condition = str(i)+'_density_'+str(j)+'_' - print(condition) - - #Division matrix - data = pd.read_excel(path+'divisions_'+str(condition)+r'.xlsx') - time = data.index.values - - #Circadian signal - data1 = pd.read_csv(path2+condition+channel1+r'_filtered.csv', index_col=0) - - imt = {str(num): [] for num in range(2, 6)} - imt_all = [] - - for col in data: - if (data1[col].dropna()).index.values[0] == 0:# and (data1[col].dropna()).index.values[-1]>210: - division_times = ((np.where(data[col]==1)[0])) - - imt_single = [] - if 2 <= len(division_times) <= 6: - for ii in range(len(division_times)-1): - imt_single.append((division_times[ii+1]-division_times[ii])*0.5) - imt_all.append((division_times[ii+1]-division_times[ii])*0.5) - imt[str(len(division_times))].append(np.mean(imt_single)) - - imt_all_dose[j] = imt_all - - # df_imt = pd.DataFrame.from_dict(imt, orient='index') - # df_imt = df_imt.transpose() - # df_imt.columns = [str(num)+' divisions' for num in range(2, 6)] - - # fig = plt.figure(figsize=(10,10)) - # df_imt.boxplot(showfliers=False) - # plt.ylabel('IMT [hours]') - # # plt.savefig(output+'IMT_boxplot_numdivisions_high_density.svg') - # plt.show() - - max_imt = np.max([max(imt_all_dose[j]) for j in dose]) - bin_width = max_imt/40 - - plt.figure(figsize=(10,8)) - for ii in dose: - print(len(imt_all_dose[ii])) - print(np.mean(imt_all_dose[ii]), np.std(imt_all_dose[ii])) - sns.histplot(imt_all_dose[ii], kde=True, stat='density', label= str(ii), color=colors[ii], bins=np.arange(16, max_imt + bin_width, bin_width)) - plt.xlabel('IMT [hours]') - plt.legend(loc='best') - # plt.savefig(output+'IMT_distributions_high_density.svg') - plt.show() - - -#%% Changing density - -density = ['high','medium','low'] - -imt_all_density = {} - -for i in density: - df = pd.DataFrame() - - condition = str(i)+'_density_'+str(dose[0])+'_' - print(condition) - - #Division matrix - data = pd.read_excel(path+'divisions_'+str(condition)+r'.xlsx') - time = data.index.values - - #Circadian signal - data1 = pd.read_csv(path2+condition+channel1+r'_filtered.csv', index_col=0) - - imt = {str(num): [] for num in range(2, 6)} - imt_all = [] - - for col in data: - if (data1[col].dropna()).index.values[0] == 0:#and (data1[col].dropna()).index.values[-1]>210: - division_times = ((np.where(data[col]==1)[0])) - - imt_single = [] - if 2 <= len(division_times) <= 6: - for ii in range(len(division_times)-1): - imt_single.append((division_times[ii+1]-division_times[ii])*0.5) - imt_all.append((division_times[ii+1]-division_times[ii])*0.5) - imt[str(len(division_times))].append(np.mean(imt_single)) - - imt_all_density[i] = imt_all -# print(imt_all_density) - -colors = {'high':'gold', 'medium':'tab:orange', 'low':'tab:blue'} - -max_imt = np.max([max(imt_all_density[i]) for i in density]) -bin_width = max_imt/40 - -plt.figure(figsize=(10,8)) -for ii in density: - print(len(imt_all_density[ii])) - print(np.mean(imt_all_density[ii]), np.std(imt_all_density[ii])) - sns.histplot(imt_all_density[ii], kde=True, stat='density', label= str(ii), color=colors[ii], bins=np.arange(16, max_imt + bin_width, bin_width)) -plt.xlabel('IMT [hours]') -plt.legend(loc='best') -# plt.savefig(output+'IMT_distributions_changing_density.svg') -plt.show() - \ No newline at end of file diff --git a/Fig4_S4/Proliferation_plots.py b/Fig4_S4/Proliferation_plots.py deleted file mode 100644 index 5c996f2..0000000 --- a/Fig4_S4/Proliferation_plots.py +++ /dev/null @@ -1,118 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:55:30 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -from mpl_toolkits.axes_grid1 import make_axes_locatable -from matplotlib.colors import LinearSegmentedColormap - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -def create_custom_colormap(num_colors): #input: maximum number divisions - all_colors = [ - (135/255, 206/255, 250/255), # skyblue - # (255/255, 228/255, 181/255), # moccasin - (255/255, 228/255, 225/255), #mistyrose - (216/255, 191/255, 216/255), # thistle - (240/255, 128/255, 128/255), # lightcoral - (95/255, 158/255, 160/255), # cadetblue - (255/255, 215/255, 0/255), # gold - ] - - if num_colors not in range(3, 7): - raise ValueError("Number of colors must be between 3 and 6 inclusive") - - colors = all_colors[:num_colors] - - cmap_name = 'custom_colormap' - return LinearSegmentedColormap.from_list(cmap_name, colors, N=256) - - -path = '/Users/nicagutu/Nextcloud/Manuscripts/Circadian-CellCycle/Data/1stExp_ilastik/Division_matrix/' -output = '/Users/nicagutu/Nextcloud/Manuscripts/Circadian-CellCycle/Figures/Fig4/' - -dose = ['untreated', '5uM', '10uM'] -density = ['high', 'medium', 'low'] -channel2 = 'cell_cycle' - -dt = 0.5 # the sampling interval, 0.5hours -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -for i in density: - for j in dose: - - condition = str(i)+'_density_'+str(j)+'_' - print(condition) - - #Division matrix - data = pd.read_excel(path+'divisions_'+str(condition)+r'.xlsx') - - IDs = data.columns - time_points = len(data.index.values) - properties = {'num_divisions':[],'time_1st':[]} - - for col in IDs: - properties['num_divisions'].append(sum(data[col])) - count=0 - for jj in range(len(data[col])): - if sum(data[col]) == 0 and count <1: - count+=1 - properties['time_1st'].append(0) - elif sum(data[col]) != 0: - if data[col][jj] == 1 and count <1: - properties['time_1st'].append(data.index.values[jj]) - count+=1 - - properties['num_divisions'],properties['time_1st'],IDs=zip(*sorted(zip(properties['num_divisions'],properties['time_1st'],IDs))) - - divisions = data.reindex(columns=IDs) - division_profile = (divisions.to_numpy()) - division_profile = division_profile.T - - for xx in range(len(IDs)): - divisions = [] - for yy in range(time_points): - if division_profile[xx, yy] == 1: - divisions.append(yy) - - if len(divisions) > 0: - for n in range(len(divisions) - 1): - division_profile[xx, divisions[n]+1:divisions[n+1]] = n+1 - - division_profile[xx, divisions[-1]+1:time_points] = len(divisions) - - if len(divisions) > 1: - division_profile[xx, divisions[-1]] = len(divisions) - - plt.figure(figsize=(10, 20)) - positions = (0, 100, 200) - label_pos = (0, 50, 100) - plt.xticks(positions, label_pos) - color_map = plt.imshow(division_profile) - color_map.set_cmap(create_custom_colormap(max(properties['num_divisions'])+1)) - plt.xlabel('Time(h)') - plt.ylabel('Cell number') - ax = plt.gca() - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.1) - cb = plt.colorbar(cax=cax) - labels_list = np.arange(0, sum(divisions), 1) - loc = labels_list+0 - cb.set_ticks(loc) - cb.set_ticklabels(labels_list) - cb.set_label('No. divisions', rotation=270, labelpad=25) - plt.savefig(output+'Division_profile_'+str(condition)+'.svg', format='svg') - plt.show() - - \ No newline at end of file diff --git a/Fig4_S4/cell_growth_plot.py b/Fig4_S4/cell_growth_plot.py deleted file mode 100644 index 7d9f06f..0000000 --- a/Fig4_S4/cell_growth_plot.py +++ /dev/null @@ -1,123 +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 = 'Input/' -output = 'Figures/Supplementary_figures/' - -# condition_combine_dict = {101:'high_density_untreated',81:'high_density_0uM',61:'high_density_1.25uM', -# 41:'high_density_2.5uM',21:'high_density_5uM',1:'high_density_10uM', -# 121:'medium_density_untreated',141:'medium_density_0uM',161:'medium_density_1.25uM', -# 181:'medium_density_2.5uM',201:'medium_density_5uM',221:'medium_density_10uM', -# 341:'low_density_untreated',321:'low_density_0uM',301:'low_density_1.25uM', -# 281:'low_density_2.5uM',261:'low_density_5uM',241:'low_density_10uM'} - -# 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',1:'high_density_10uM'} - -# 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) - -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