############################################################################### # Info about how the model generally works: # 1) Equilibration: Desired Mcl-1 concentration is distributed between cytoplasm and mitochondria # of desired volume (This is not part of the simulation results) # 2) G1: G1 phase is simulated with its respective parameters, Concentrations and volumes are saved # and passed on to the next phase # 3) S/G2: S/G2 phase are simulated as one phase with respective parameters. The end values from G1 # are taken as input and are simulated further, saved, and passed on to he next phase # 4) M: M phase is simulated from end point of S/G2 with respective parameters. End values are saved.5 # 5) G1 (second time): Concentrations from end of M phase are kept, but volume is halved (thereby also # protein number, while concentrations retains. # 6) Cells cycle for as many cycles as defined in input parameters. # 7) After the simulation, array of concentrations and volumes over time are saved for each cell. ############################################################################### # Info about cell to cell heterogeneity implementation # 1) Heterogeneity in Mcl-1 distribution through random sampling from experimentally observed Mcl-1 # distribution heterogeneity # 2) Heterogeneity in Mcl-1 expression through random sampling from experimentally observed Mcl-1 # expression heterogeneity1 # 3) Heterogeneity in cell size through random sampling from experimentally observed cell sizes # 4) Heterogeneity in cell cycle stage through choosing random starting points for the simulation # --> All cells are simulated for one more cycle as defined in the input parameters # --> A random timepoint of one cell cycle duration is selected for each cell # --> The arrays over time are cut, so that t=0 equals that starting point and also end point # after n cycles ################################################################################ # Info about code structure # 1) Cell cycle phase models are defined first # 2) All functions for data handling, saving simulations and passing values to the next simulation etc. # are defined # 3) Parameters of the simulation are defined # 4) Loop structure to simulate cells, process data and save results, optionally also in batches to not # exceed memory issues with larger simulations # Import packages import tellurium as te import roadrunner import matplotlib.pyplot as plt import matplotlib from matplotlib.gridspec import GridSpec import numpy as np import xlsxwriter import csv import math import os from datetime import datetime import seaborn as sns from sklearn.metrics import auc import pandas as pd import random import time import shutil import ipykernel import requests from IPython.display import display, Javascript from notebook.notebookapp import list_running_servers # Load eq model # This part creates the virtual cell with cytoplasmic and mitochondrial compartment # according to the experimental input volumes # It takes the initial Mcl-1 concentration and distributes it according # to the experimental k_dist between cytoplasm and mitochondria def load_eq(): ModelName = "Model for equilibration" global toymodel_eq toymodel_eq = (''' model toymodel_eq() //------------------------------------------------------------------------ // Input /* Mito_ini and Cyto_ini are needed, since variable compartments like Mito, Cyto or Cell can't be assigned with an initial value from outside of the model definition. Therefore, Mito, Cyto and Cell volumes are defined based on the initial values of Mito_ini and Cyto_ini */ compartment Mito_ini = 0; # µm^3 compartment Cyto_ini = 0; # µm^3 var compartment Mito in Cyto := Mito_ini; # µm^3 var compartment Cyto in Cell := Cyto_ini; # µm^3 var compartment Cell := Mito+Cyto; # µm^3 var substanceOnly species Mcl1 in Cell = 0; # nM var substanceOnly species Mcl1M in Mito = 0; # nM var substanceOnly species Mcl1C in Cyto = 0; # nM //------------------------------------------------------------------------ // Reactions Mcl1 => Mcl1C; distr*Mcl1; # Distribution of initial Mcl1 in Cell Mcl1C => Mcl1M; k_on*(Mcl1C/Cyto); # Mcl1 from Cyto onto Mito Mcl1M => Mcl1C; k_off*(Mcl1M/Mito); # Mcl1 from Mito into Cyto //------------------------------------------------------------------------ //Parameters distr = 0.0005; # s^-1 k_d = 0; k_off = 0; # s^-1 k_on := k_off/k_d; # s^-1 end ''') toymodel_eq = te.loada(toymodel_eq) ODEs_eq = te.getODEsFromModel(toymodel_eq) #toymodel_eq.exportToSBML("toymodel_eq.sbml") # Load Sim_G1 model # This part models the G1 phase of the cell cycle. # Mitochondrial and cytoplasmic compartments increase in size according to growth rates from literature # Mcl-1 is degraded according to its half-life (mass-action) and produced in a similar rate # Mcl-1 production increases exponentially, similar to compartment growth # Degradation/Synthesis as well as k_on and k_off values can be adjusted from outside the function # k_on is chosen artifically (for now) # k_off decreases exponentially with progression in cell cycle to mimic mitochondrial accumulation of Mcl-1 def load_g1(): ModelName = "Model for simulating G1" global toymodel_sim_g1 toymodel_sim_g1 = (''' model toymodel_sim_g1() //------------------------------------------------------------------------ // Input /* Mito_ini and Cyto_ini are needed, since variable compartments like Mito, Cyto or Cell can't be assigned with an initial value from outside of the model definition. Therefore, Mito, Cyto and Cell volumes are defined based on the initial values of Mito_ini and Cyto_ini */ compartment Mito_ini = 0; # µm^3 compartment Cyto_ini = 0; # µm^3 var compartment Mito in Cyto := Mito_ini*exp(growthrate*(time)); # µm^3 var compartment Cyto in Cell := Cyto_ini*exp(growthrate*(time)); # µm^3 var compartment Cell := Mito+Cyto; # µm^3 var substanceOnly species Mcl1M in Mito = 0; # nM var substanceOnly species Mcl1C in Cyto = 0; # nM //------------------------------------------------------------------------ // Reactions Mcl1C => Mcl1M; k_on*(Mcl1C/Cyto); # Mcl1 from Cyto onto Mito Mcl1M => Mcl1C; k_off*(Mcl1M/Mito); # Mcl1 from Mito into Cyto Mcl1C =>; k_deg*Mcl1C*f_deg; # Mcl1 Degradation in Cyto => Mcl1C; k_syn*syn_test*f_syn; # Mcl1 Synthesis in Cyto syn_test := k_mcl1eq*exp(growthrate*(time)) # to debug and access synthesis factor deg_test := k_deg*Mcl1C*f_deg # to access degradation rate //------------------------------------------------------------------------ //Parameters kd = 0; kd_test := kd*(2-(0.35/(1+exp(-0.0009*(time-3500)))-1.7*10^(-15)*(time+10000)^(3)+0.99)); #kd_test := kd k_off = 0; k_on := (k_off/kd)*(2-(0.35/(1+exp(-0.0009*(time-3500)))-1.7*10^(-15)*(time+10000)^(3)+0.99)); #k_on := k_off/kd kd_rate = 0 k_deg = 0; k_syn = 0; f_deg = 0; f_syn = 0; k_mcl1eq = 0; t_treatment = 0; chx = 0; chx_effect = 0; btz = 0; btz_effect = 0; growthrate = 0; //------------------------------------------------------------------------ //Events at(((time)>=t_treatment)&&(chx==1)): k_syn = k_syn*chx_effect; # stops Mcl-1 synthesis, if desired at(((time)>=t_treatment)&&(btz==1)): k_deg = k_deg*btz_effect; # stops Mcl-1 degradation, if desired end ''') toymodel_sim_g1 = te.loada(toymodel_sim_g1) ODEs_sim_g1 = te.getODEsFromModel(toymodel_sim_g1) #print("ODEs Toymodel_sim_g1") #print(ODEs_sim_g1) #toymodel_sim_g1.exportToSBML("toymodel_sim_g1.sbml") # Load Sim_S/G2 model # This part models the S/G2 phase of the cell cycle. # Mitochondrial and cytoplasmic compartments increase in size according to growth rates from literature # Mcl-1 is degraded according to its half-life (mass-action) and produced in a similar rate # Mcl-1 production increases exponentially, similar to compartment growth # Degradation/Synthesis as well as k_on and k_off values can be adjusted from outside the function # k_on is chosen artifically (for now) # k_off decreases exponentially with progression in cell cycle to mimic mitochondrial accumulation of Mcl-1 def load_sg2(): ModelName = "Model for simulating S/G2" global toymodel_sim_sg2 toymodel_sim_sg2 = (''' model toymodel_sim_sg2() //------------------------------------------------------------------------ // Input /* Mito_ini and Cyto_ini are needed, since variable compartments like Mito, Cyto or Cell can't be assigned with an initial value from outside of the model definition. Therefore, Mito, Cyto and Cell volumes are defined based on the initial values of Mito_ini and Cyto_ini */ compartment Mito_ini = 0; # µm^3 compartment Cyto_ini = 0; # µm^3 var compartment Mito in Cyto := Mito_ini*exp(growthrate*(time)); # µm^3 var compartment Cyto in Cell := Cyto_ini*exp(growthrate*(time)); # µm^3 var compartment Cell := Mito+Cyto; # µm^3 var substanceOnly species Mcl1M in Mito = 0; # nM var substanceOnly species Mcl1C in Cyto = 0; # nM //------------------------------------------------------------------------ // Reactions Mcl1C => Mcl1M; k_on*(Mcl1C/Cyto); # Mcl1 from Cyto onto Mito Mcl1M => Mcl1C; k_off*(Mcl1M/Mito); # Mcl1 from Mito into Cyto Mcl1C =>; k_deg*Mcl1C*f_deg*(2-exp(0.000013*((time)))); # Mcl1 Degradation in Cyto #Mcl1C =>; k_deg*Mcl1C*f_deg; # Mcl1 Degradation in Cyto => Mcl1C; k_syn*syn_test*f_syn; # Mcl1 Synthesis in Cyto syn_test := k_mcl1eq*exp(growthrate*((time))) # to debug and save synthesis factor deg_test := k_deg*Mcl1C*f_deg # to get live access to degradation rate //------------------------------------------------------------------------ //Parameters kd = 0; kd_test := kd*(1.578-(0.1/(1+exp(0.0005*(time-3500)))-1.918859*10^(-15)*(time+10000)^(3)+0.906)); #kd_test := kd k_off = 0; k_on := k_off/kd*(1.578-(0.1/(1+exp(0.0005*(time-3500)))-1.918859*10^(-15)*(time+10000)^(3)+0.906)); #k_on := k_off/kd kd_rate = 0 k_deg = 0; k_syn = 0; f_deg = 0; f_syn = 0; k_mcl1eq = 0; t_treatment = 0; chx = 0; chx_effect = 0; btz = 0; btz_effect = 0; growthrate = 0; //------------------------------------------------------------------------ //Events at(((time)>=t_treatment)&&(chx==1)): k_syn = k_syn*chx_effect; at(((time)>=t_treatment)&&(btz==1)): k_deg = k_deg*btz_effect; end ''') toymodel_sim_sg2 = te.loada(toymodel_sim_sg2) ODEs_sim_sg2 = te.getODEsFromModel(toymodel_sim_sg2) #print("ODEs Toymodel_sim_sg2") #print(ODEs_sim_sg2) #toymodel_sim_sg2.exportToSBML("toymodel_sim_sg2.sbml") # Load M model # This part models the G1 phase of the cell cycle. # Mitochondrial and cytoplasmic compartments do not increase in size anymore # Mcl-1 is degraded according to its half-life (mass-action) and produced in a similar rate # Mcl-1 production increases exponentially, similar to compartment growth # Degradation/Synthesis as well as k_on and k_off values can be adjusted from outside the function # k_on is chosen artifically (for now) # k_off decreases exponentially with progression in cell cycle to mimic mitochondrial accumulation of Mcl-1 def load_m(): ModelName = "Model for simulating M" global toymodel_sim_m toymodel_sim_m = (''' model toymodel_sim_m() //------------------------------------------------------------------------ // Input /* Mito_ini and Cyto_ini are needed, since variable compartments like Mito, Cyto or Cell can't be assigned with an initial value from outside of the model definition. Therefore, Mito, Cyto and Cell volumes are defined based on the initial values of Mito_ini and Cyto_ini */ compartment Mito_ini = 0; # µm^3 compartment Cyto_ini = 0; # µm^3 var compartment Mito in Cyto := Mito_ini*exp(growthrate*(time));; # µm^3 var compartment Cyto in Cell := Cyto_ini*exp(growthrate*(time));; # µm^3 var compartment Cell := Mito+Cyto; # µm^3 var substanceOnly species Mcl1M in Mito = 0; # nM var substanceOnly species Mcl1C in Cyto = 0; # nM //------------------------------------------------------------------------ // Reactions Mcl1C => Mcl1M; k_on*(Mcl1C/Cyto); # Mcl1 from Cyto onto Mito Mcl1M => Mcl1C; k_off*(Mcl1M/Mito); # Mcl1 from Mito into Cyto Mcl1C =>; k_deg*Mcl1C*f_deg; # Mcl1 Degradation in Cyto => Mcl1C; k_syn*syn_test*f_syn; # Mcl1 Synthesis in Cyto syn_test := k_mcl1eq*exp(growthrate*((time))) deg_test := k_deg*Mcl1C*f_deg //------------------------------------------------------------------------ //Parameters kd = 0; # 1.58 and + instead of mines kd_test := kd*(1.75-(0.14/(1+exp(0.0005228942*(time-33500-33048)))-1.8006859*10^(-15)*(time+10000-33048)^(3)+0.906)); #kd_test := kd k_off = 0; k_on := k_off/kd*(1.75-(0.14/(1+exp(0.0005228942*(time-33500-33048)))-1.8006859*10^(-15)*(time+10000-33048)^(3)+0.906)); #k_on := k_off/kd kd_rate = 0; k_deg = 0; k_syn = 0; f_deg = 0; f_syn = 0; k_mcl1eq = 0; t_treatment = 0; chx = 0; chx_effect = 0; btz = 0; btz_effect = 0; growthrate = 0; //------------------------------------------------------------------------ //Events at(((time)>=t_treatment)&&(chx==1)): k_syn = k_syn*chx_effect; at(((time)>=t_treatment)&&(btz==1)): k_deg = k_deg*btz_effect; end ''') toymodel_sim_m = te.loada(toymodel_sim_m) ODEs_sim_m = te.getODEsFromModel(toymodel_sim_m) #print("ODEs Toymodel_sim_m") #print(ODEs_sim_m) #toymodel_sim_m.exportToSBML("toymodel_sim_m.sbml") # eq_preparation broadcasts defined input values into the antimony model def eq_preparation(): # Array of all predefined input values Input_EQ = [] Input_EQ.append(float(v_mito_ini)*0.6022) # Mito_vol in [µm^3], avogardo-corrected by 0.6022, v_mito_ini is defined in loop structure Input_EQ.append(float(v_cyto_ini)*0.6022) # Cyto_vol in [µm^3], avogardo-corrected by 0.6022, v_cyto_ini is defined in loop structure Input_EQ.append(float(c_mcl1_ini)) # Mcl1 in the cell [nM], c_mcl1_ini is defined in loop structure Input_EQ.append(float(kd_eq)) # k_d cyto/mito for equilibration Input_EQ.append(float(k_off_eq)) # k_on (arbitrarely) of Mcl-1 on mitochondria # Values are defined within the loaded model toymodel_eq.reset() toymodel_eq.Mito_ini = (float(Input_EQ[0])) # Set initial, total Mcl1 concentration toymodel_eq.Cyto_ini = (float(Input_EQ[1])) # Set initial, total Mcl1 concentration toymodel_eq.Mcl1 = (float(Input_EQ[2])) # Set initial, total Mcl1 concentration toymodel_eq.k_d = (float(Input_EQ[3])) # set initial k_d toymodel_eq.k_off = (float(Input_EQ[4])) # set initial k_on # Input values for G1 phase def g1_preparation(t_time, treat, time_passed, treat_con, treat_dif): temp_time = t_time - treat_dif # defines temporary time variable - at which time in the sub-simulation the treatment should start if treat_con == 1: # if treatment has already started, temp_time = 1 # treatemnt starts from the beginning of the sub-simulation Input_G1 = [] Input_G1.append(float(vol_mito_g1)*0.6022) # [0] = Equilibrated Mito volume, saved as last value from EQ Input_G1.append(float(vol_cyto_g1)*0.6022) # [1] = Equilibrated Cyto volume, saved as last value from EQ Input_G1.append(float(Mcl1M_g1)) # [2] = Equilibrated Mito Mcl1 concentration, saved as last value from EQ Input_G1.append(float(Mcl1C_g1)) # [3] = Equilibrated Cyto Mcl1 concentration, saved as last value from EQ Input_G1.append(float(Mcl1C_part_eq_g1)) # [4] = Equilibrated Mcl1 amount for Synthesis rate, saved as last value from EQ Input_G1.append(float(mcl1_deg_g1)) # [5] = k_deg rate Input_G1.append(float(mcl1_syn_g1)) # [6] = k_syn rate Input_G1.append(float(temp_time)) # [7] = Treatment time Input_G1.append(float(chx*treat)) # [8] = CHX treatment? 0 = no, 1 = yes Input_G1.append(float(btz*treat)) # [9] = BTZ treatment? 0 = no, 1 = yes Input_G1.append(float(chx_effect)) # [10] = CHX effect from 0 (100 %) to 1 (0 %) Input_G1.append(float(btz_effect)) # [11] = BTZ effect from 0 (100 %) to 1 (0 %) Input_G1.append(float(growth_g1)) # [12] = Cell growth rate in G1 (1/V in s^-1) Input_G1.append(float(kd_g1)) # [13] = KD Mito/Cyto in G1 Input_G1.append(float(kd_rate_g1)) # [14] = KD rate Mito/Cyto in G1 Input_G1.append(float(deg_g1)) # [15] = degradation factor Input_G1.append(float(syn_g1)) # [16] = synthesis factor Input_G1.append(float(k_off)) # [17] = k_on toymodel_sim_g1.reset() toymodel_sim_g1.Mito_ini = (float(Input_G1[0])) # Set initial, total Mcl1 concentration toymodel_sim_g1.Cyto_ini = (float(Input_G1[1])) # Set initial, total Mcl1 concentration toymodel_sim_g1.Mcl1M = (float(Input_G1[2])) # Set equilibrated Mcl1M toymodel_sim_g1.Mcl1C = (float(Input_G1[3])) # Set equilibrated Mcl1C toymodel_sim_g1.k_mcl1eq = (float(Input_G1[4])) # Set equilibrated Mcl1C amount for synthesis toymodel_sim_g1.k_deg = (float(Input_G1[5])) # Set k_deg toymodel_sim_g1.k_syn = (float(Input_G1[6])) # Set k_syn toymodel_sim_g1.t_treatment = (float(Input_G1[7])) # Set treatment time toymodel_sim_g1.chx = (float(Input_G1[8])) # CHX treatment? toymodel_sim_g1.btz = (float(Input_G1[9])) # BTZ treatment? toymodel_sim_g1.chx_effect = (float(Input_G1[10])) # CHX effect toymodel_sim_g1.btz_effect = (float(Input_G1[11])) # BTZ effect toymodel_sim_g1.growthrate = (float(Input_G1[12])) # Cell growth toymodel_sim_g1.kd = (float(Input_G1[13])) # KD Mcl1 toymodel_sim_g1.kd_rate = (float(Input_G1[14])) # KD rate Mcl1 toymodel_sim_g1.f_deg = (float(Input_G1[15])) # deg factor toymodel_sim_g1.f_syn = (float(Input_G1[16])) # syn factor toymodel_sim_g1.k_off = (float(Input_G1[17])) # k_on # Input values for S/G2 phase def sg2_preparation(t_time, treat, time_passed, treat_con, treat_dif): temp_time = t_time - treat_dif if treat_con == 1: temp_time = 1 Input_SG2 = [] Input_SG2.append(float(vol_mito_sg2)*0.6022) # [0] = Equilibrated Mito volume Input_SG2.append(float(vol_cyto_sg2)*0.6022) # [1] = Equilibrated Cyto volume Input_SG2.append(float(Mcl1M_sg2)) # [2] = Equilibrated Mito Mcl1 concentration Input_SG2.append(float(Mcl1C_sg2)) # [3] = Equilibrated Cyto Mcl1 concentration Input_SG2.append(float(Mcl1C_part_eq_sg2)) # [4] = Equilibrated Mcl1 amount for Synthesis rate Input_SG2.append(float(mcl1_deg_sg2)) # [5] = k_deg rate/2 Input_SG2.append(float(mcl1_syn_sg2)) # [6] = k_syn rate Input_SG2.append(float(temp_time)) # [7] = Treatment time Input_SG2.append(float(chx*treat)) # [8] = CHX treatment? 0 = no, 1 = yes Input_SG2.append(float(btz*treat)) # [9] = BTZ treatment? 0 = no, 1 = yes Input_SG2.append(float(chx_effect)) # [10] = CHX effect from 0 (100 %) to 1 (0 %) Input_SG2.append(float(btz_effect)) # [11] = BTZ effect from 0 (100 %) to 1 (0 %) Input_SG2.append(float(growth_sg2)) # [12] = Cell growth rate in G1 (1/V in s^-1) Input_SG2.append(float(kd_sg2)) # [13] = KD Mito/Cyto in SG2 Input_SG2.append(float(kd_rate_sg2)) # [14] = KD rate Mito/Cyto in SG2 Input_SG2.append(float(deg_sg2)) # [15] = degradation factor Input_SG2.append(float(syn_sg2)) # [16] = synthesis factor Input_SG2.append(float(k_off)) # [17] = k_on toymodel_sim_sg2.reset() toymodel_sim_sg2.Mito_ini = (float(Input_SG2[0])) # Set initial, total Mcl1 concentration toymodel_sim_sg2.Cyto_ini = (float(Input_SG2[1])) # Set initial, total Mcl1 concentration toymodel_sim_sg2.Mcl1M = (float(Input_SG2[2])) # Set equilibrated Mcl1M toymodel_sim_sg2.Mcl1C = (float(Input_SG2[3])) # Set equilibrated Mcl1C toymodel_sim_sg2.k_mcl1eq = (float(Input_SG2[4])) # Set equilibrated Mcl1C amount for synthesis toymodel_sim_sg2.k_deg = (float(Input_SG2[5])) # Set k_deg toymodel_sim_sg2.k_syn = (float(Input_SG2[6])) # Set k_syn toymodel_sim_sg2.t_treatment = (float(Input_SG2[7])) # Set treatment time toymodel_sim_sg2.chx = (float(Input_SG2[8])) # CHX treatment? toymodel_sim_sg2.btz = (float(Input_SG2[9])) # BTZ treatment? toymodel_sim_sg2.chx_effect = (float(Input_SG2[10])) # CHX effect toymodel_sim_sg2.btz_effect = (float(Input_SG2[11])) # BTZ effect toymodel_sim_sg2.growthrate = (float(Input_SG2[12])) # Cell growth toymodel_sim_sg2.kd = (float(Input_SG2[13])) # KD Mcl1 toymodel_sim_sg2.kd_rate = (float(Input_SG2[14])) # KD rate Mcl1 toymodel_sim_sg2.f_deg = (float(Input_SG2[15])) # deg factor toymodel_sim_sg2.f_syn = (float(Input_SG2[16])) # syn factor toymodel_sim_sg2.k_off = (float(Input_SG2[17])) # k_on # Input values for M phase def m_preparation(t_time, treat, time_passed, treat_con, treat_dif): temp_time = t_time - treat_dif if treat_con == 1: temp_time = 1 Input_M = [] Input_M.append(float(vol_mito_m)*0.6022) # [0] = Equilibrated Mito volume Input_M.append(float(vol_cyto_m)*0.6022) # [1] = Equilibrated Cyto volume Input_M.append(float(Mcl1M_m)) # [2] = Equilibrated Mito Mcl1 concentration Input_M.append(float(Mcl1C_m)) # [3] = Equilibrated Cyto Mcl1 concentration Input_M.append(float(Mcl1C_part_eq_m)) # [4] = Equilibrated Mcl1 amount for Synthesis rate Input_M.append(float(mcl1_deg_m)) # [5] = k_deg rate/3 Input_M.append(float(mcl1_syn_m)) # [6] = k_syn rate Input_M.append(float(temp_time)) # [7] = Treatment time Input_M.append(float(chx*treat)) # [8] = CHX treatment? 0 = no, 1 = yes Input_M.append(float(btz*treat)) # [9] = BTZ treatment? 0 = no, 1 = yes Input_M.append(float(chx_effect)) # [10] = CHX effect from 0 (100 %) to 1 (0 %) Input_M.append(float(btz_effect)) # [11] = BTZ effect from 0 (100 %) to 1 (0 %) Input_M.append(float(growth_m)) # [12] = Cell growth rate in G1 (1/V in s^-1) Input_M.append(float(kd_m)) # [13] = KD Mito/Cyto in M Input_M.append(float(kd_rate_m)) # [14] = KD rate Mito/Cyto in M Input_M.append(float(deg_m)) # [15] = degradation factor Input_M.append(float(syn_m)) # [16] = synthesis factor Input_M.append(float(k_off)) # [17] = k_on toymodel_sim_m.reset() toymodel_sim_m.Mito_ini = (float(Input_M[0])) # Set initial, total Mcl1 concentration toymodel_sim_m.Cyto_ini = (float(Input_M[1])) # Set initial, total Mcl1 concentration toymodel_sim_m.Mcl1M = (float(Input_M[2])) # Set equilibrated Mcl1M toymodel_sim_m.Mcl1C = (float(Input_M[3])) # Set equilibrated Mcl1C toymodel_sim_m.k_mcl1eq = (float(Input_M[4])) # Set equilibrated Mcl1C amount for synthesis toymodel_sim_m.k_deg = (float(Input_M[5])) # Set k_deg toymodel_sim_m.k_syn = (float(Input_M[6])) # Set k_syn toymodel_sim_m.t_treatment = (float(Input_M[7])) # Set treatment time toymodel_sim_m.chx = (float(Input_M[8])) # CHX treatment? toymodel_sim_m.btz = (float(Input_M[9])) # BTZ treatment? toymodel_sim_m.chx_effect = (float(Input_M[10])) # CHX effect toymodel_sim_m.btz_effect = (float(Input_M[11])) # BTZ effect toymodel_sim_m.growthrate = (float(Input_M[12])) # Cell growth toymodel_sim_m.kd = (float(Input_M[13])) # KD Mcl1 toymodel_sim_m.kd_rate = (float(Input_M[14])) # KD rate Mcl1 toymodel_sim_m.f_deg = (float(Input_M[15])) # deg factor toymodel_sim_m.f_syn = (float(Input_M[16])) # syn factor toymodel_sim_m.k_off = (float(Input_M[17])) # k_on def plot_results_eq(n): # Plot Results #toymodel_eq.plot(xtitle="time [s]", ytitle="concentration [nM]", title="Mcl-1 Equilibration") #te.plotArray(sim1) plt.figure() plt.plot(sim_eq[:,0],sim_eq[:,1]/0.6022, label = "Mito Volume") plt.plot(sim_eq[:,0],sim_eq[:,2]/0.6022, label = "Cyto Volume") plt.plot(sim_eq[:,0],sim_eq[:,3]/0.6022, label = "Cell Volume") plt.xlabel("time [s]") plt.ylabel("Volume [µm^3]") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Compartment volumes over time of Cell Nr. " +str(n)) plt.show() plt.figure() plt.plot(sim_eq[:,0],sim_eq[:,4], label = "Initital Mcl1") plt.plot(sim_eq[:,0],sim_eq[:,5], label = "Mitochondrial Mcl1") plt.plot(sim_eq[:,0],sim_eq[:,6], label = "Cytoplasmic Mcl1") plt.xlabel("time [s]") plt.ylabel("Concentration [nM]") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Species concentrations over time of Cell Nr. " +str(n)) plt.show() plt.figure() plt.plot(sim_eq[:,0],sim_eq[:,7], label = "Initital Mcl1 particles") plt.plot(sim_eq[:,0],sim_eq[:,8], label = "Mitochondrial Mcl1 particles") plt.plot(sim_eq[:,0],sim_eq[:,9], label = "Cytoplasmic Mcl1 particles") plt.xlabel("time [s]") plt.ylabel("Particles") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Species particles over time of Cell Nr. " +str(n)) plt.show def save_results_eq(n): global vol_mito_eq global vol_cyto_eq global vol_cell_eq vol_mito_eq = sim_eq[Time_Simulation_EQ[2],1]/0.6022 # mito volume at last point of equilibration vol_cyto_eq = sim_eq[Time_Simulation_EQ[2],2]/0.6022 # cyto volume at last point of equilibration vol_cell_eq = sim_eq[Time_Simulation_EQ[2],3]/0.6022 # cell volume at last point of equilibration global Mcl1I_eq global Mcl1M_eq global Mcl1C_eq Mcl1I_eq = sim_eq[Time_Simulation_EQ[2],4] # Initial Mcl-1 conc. at last point of equilibration Mcl1M_eq = sim_eq[Time_Simulation_EQ[2],5] # Mito Mcl-1 conc. at last point of equilibration Mcl1C_eq = sim_eq[Time_Simulation_EQ[2],6] # Cyto Mcl-1 conc. at last point of equilibration global Mcl1I_part_eq_g1 global Mcl1M_part_eq_g1 global Mcl1C_part_eq_g1 Mcl1I_part_eq_g1 = sim_eq[Time_Simulation_EQ[2],7] # Initial Mcl-1 particles at last point of equilibration Mcl1M_part_eq_g1 = sim_eq[Time_Simulation_EQ[2],8] # Mito Mcl-1 particles at last point of equilibration Mcl1C_part_eq_g1 = sim_eq[Time_Simulation_EQ[2],9] # Cyto Mcl-1 particles at last point of equilibration # Grab equilibrated volumes and concentration from Toymodel_EQ # Initial values for SG2 and M are not needed, variables must be defined # --> They are changed later on global vol_mito_g1 global vol_cyto_g1 global vol_mito_sg2 global vol_cyto_sg2 global vol_mito_m global vol_cyto_m vol_mito_g1 = vol_mito_eq vol_cyto_g1 = vol_cyto_eq vol_mito_sg2 = vol_mito_eq vol_cyto_sg2 = vol_cyto_eq vol_mito_m = vol_mito_eq vol_cyto_m = vol_cyto_eq global Mcl1M_g1 global Mcl1C_g1 global Mcl1M_sg2 global Mcl1C_sg2 global Mcl1M_m global Mcl1C_m Mcl1M_g1 = Mcl1M_eq Mcl1C_g1 = Mcl1C_eq Mcl1M_sg2 = Mcl1M_eq Mcl1C_sg2 = Mcl1C_eq Mcl1M_m = Mcl1M_eq Mcl1C_m = Mcl1C_eq temp_eq_time=(sim_eq[:,0]) temp_eq_MitoV=(sim_eq[:,1]/0.6022) temp_eq_CytoV=(sim_eq[:,2]/0.6022) temp_eq_CellV=(sim_eq[:,3]/0.6022) temp_eq_Mcl1_C=(sim_eq[:,4]) temp_eq_Mcl1M_C=(sim_eq[:,5]) temp_eq_Mcl1C_C=(sim_eq[:,6]) temp_eq_Mcl1_P=(sim_eq[:,7]) temp_eq_Mcl1M_P=(sim_eq[:,8]) temp_eq_Mcl1C_P=(sim_eq[:,9]) for t in range(len(temp_eq_time)): if t % 1000 == 0: eq_time[n].append(temp_eq_time[t]) eq_MitoV[n].append(temp_eq_MitoV[t]) eq_CytoV[n].append(temp_eq_CytoV[t]) eq_CellV[n].append(temp_eq_CellV[t]) eq_Mcl1_C[n].append(temp_eq_Mcl1_C[t]) eq_Mcl1C_C[n].append(temp_eq_Mcl1C_C[t]) eq_Mcl1M_C[n].append(temp_eq_Mcl1M_C[t]) eq_Mcl1_P[n].append(temp_eq_Mcl1_P[t]) eq_Mcl1C_P[n].append(temp_eq_Mcl1C_P[t]) eq_Mcl1M_P[n].append(temp_eq_Mcl1M_P[t]) # Save simulation results from G1 def save_results_g1(c,ccd,n): temp_time=(sim_g1[:,0]) # save time duration temp_MitoV=(sim_g1[:,1]/0.6022) # save mito volume over time temp_CytoV=(sim_g1[:,2]/0.6022) # save cyto volume over time temp_CellV=(sim_g1[:,3]/0.6022) # save cell volume over time temp_Mcl1M_C=(sim_g1[:,4]) # save Mcl-1 mito concentration over time temp_Mcl1C_C=(sim_g1[:,5]) # save Mcl-1 cyto concentration over time temp_Mcl1M_P=(sim_g1[:,6]) # save Mcl-1 mito particles over time temp_Mcl1C_P=(sim_g1[:,7]) # save Mcl-1 cyto particles over time temp_kd_new=(sim_g1[:,25]) # save kd_new values --> decreasing exponentially over time # k_off := (Cyto/Mito)*k_on*kd*(1/exp(kd_rate*(time))); # s^-1, increases over time like compartment volumes # kd_new := kd*(1/exp(kd_rate*(time))) # saves value to use in next cell cycle phase temp_syn_test=(sim_g1[:,23]) # save syn_test --> save synthesis-factor temp_deg_test=(sim_g1[:,24]) # save degradation-factor #temp_time=(temp_time)#+(Sim_start[n]) # add Sim_start time to every entry of the time array (shift time array) for t in range(len(temp_time)): if t % 1000 == 0: # take every 1000th element of every result array and save them in new result array result_time[n].append(temp_time[t]+c*ccd) # add time to final result array (+ passed cell cycles) result_MitoV[n].append(temp_MitoV[t]) result_CytoV[n].append(temp_CytoV[t]) result_CellV[n].append(temp_CellV[t]) result_Mcl1M_C[n].append(temp_Mcl1M_C[t]) result_Mcl1C_C[n].append(temp_Mcl1C_C[t]) result_Mcl1M_P[n].append(temp_Mcl1M_P[t]) result_Mcl1C_P[n].append(temp_Mcl1C_P[t]) result_kd_new[n].append(temp_kd_new[t]) result_syn_test[n].append(temp_syn_test[t]) result_deg_test[n].append(temp_deg_test[t]) global vol_mito_sg2 global vol_cyto_sg2 global Mcl1M_sg2 global Mcl1C_sg2 global Mcl1C_part_eq_sg2 global kd_sg2 vol_mito_sg2 = result_MitoV[n][-1] # starting mito volume for S/G2 phase is last value from G1 phase vol_cyto_sg2 = result_CytoV[n][-1] # starting cyto volume for S/G2 phase is last value from G1 phase Mcl1M_sg2 = result_Mcl1M_C[n][-1] # starting mito Mcl-1 concentration for S/G2 phase is last value from G1 phase Mcl1C_sg2 = result_Mcl1C_C[n][-1] # starting cyto Mcl-1 concentration for S/G2 phase is last value from G1 phase Mcl1C_part_eq_sg2 = result_syn_test[n][-1] # starting SG/2 synthesis factor is last value from G1 phase kd_sg2 = result_kd_new[n][-1] # starting SG/2 k_distr is last value from G1 phase # Save simulation results from S/G2 def save_results_sg2(c,ccd,cd1,n): temp_time=(sim_sg2[:,0]) temp_MitoV=(sim_sg2[:,1]/0.6022) temp_CytoV=(sim_sg2[:,2]/0.6022) temp_CellV=(sim_sg2[:,3]/0.6022) temp_Mcl1M_C=(sim_sg2[:,4]) temp_Mcl1C_C=(sim_sg2[:,5]) temp_Mcl1M_P=(sim_sg2[:,6]) temp_Mcl1C_P=(sim_sg2[:,7]) temp_kd_new=(sim_sg2[:,25]) temp_syn_test=(sim_sg2[:,23]) temp_deg_test=(sim_sg2[:,24]) #temp_time=(temp_time)#+(Sim_start[n]) for t in range(len(temp_time)): if t % 1000 == 0: result_time[n].append(temp_time[t]+c*ccd+cd1) result_MitoV[n].append(temp_MitoV[t]) result_CytoV[n].append(temp_CytoV[t]) result_CellV[n].append(temp_CellV[t]) result_Mcl1M_C[n].append(temp_Mcl1M_C[t]) result_Mcl1C_C[n].append(temp_Mcl1C_C[t]) result_Mcl1M_P[n].append(temp_Mcl1M_P[t]) result_Mcl1C_P[n].append(temp_Mcl1C_P[t]) result_kd_new[n].append(temp_kd_new[t]) result_syn_test[n].append(temp_syn_test[t]) result_deg_test[n].append(temp_deg_test[t]) global vol_mito_m global vol_cyto_m global Mcl1M_m global Mcl1C_m global Mcl1C_part_eq_m global kd_m vol_mito_m = result_MitoV[n][-1] vol_cyto_m = result_CytoV[n][-1] Mcl1M_m = result_Mcl1M_C[n][-1] Mcl1C_m = result_Mcl1C_C[n][-1] Mcl1C_part_eq_m = result_syn_test[n][-1] kd_m = result_kd_new[n][-1] # Save simulation results from M def save_results_m(c,ccd,cd2,n): temp_time=(sim_m[:,0]) temp_MitoV=(sim_m[:,1]/0.6022) temp_CytoV=(sim_m[:,2]/0.6022) temp_CellV=(sim_m[:,3]/0.6022) temp_Mcl1M_C=(sim_m[:,4]) temp_Mcl1C_C=(sim_m[:,5]) temp_Mcl1M_P=(sim_m[:,6]) temp_Mcl1C_P=(sim_m[:,7]) temp_syn_test=(sim_m[:,23]) temp_deg_test=(sim_m[:,24]) temp_kd_new=(sim_m[:,25]) #temp_time=(temp_time)#+(Sim_start[n]) for t in range(len(temp_time)): if t % 1000 == 0: result_time[n].append(temp_time[t]+c*ccd+cd2) result_MitoV[n].append(temp_MitoV[t]) result_CytoV[n].append(temp_CytoV[t]) result_CellV[n].append(temp_CellV[t]) result_Mcl1M_C[n].append(temp_Mcl1M_C[t]) result_Mcl1C_C[n].append(temp_Mcl1C_C[t]) result_Mcl1M_P[n].append(temp_Mcl1M_P[t]) result_Mcl1C_P[n].append(temp_Mcl1C_P[t]) result_kd_new[n].append(temp_kd_new[t]) result_syn_test[n].append(temp_syn_test[t]) result_deg_test[n].append(temp_deg_test[t]) global vol_mito_g1 global vol_cyto_g1 global Mcl1M_g1 global Mcl1C_g1 global Mcl1C_part_eq_g1 vol_mito_g1 = vol_mito_eq vol_cyto_g1 = vol_cyto_eq Mcl1M_g1 = result_Mcl1M_C[n][-1] Mcl1C_g1 = result_Mcl1C_C[n][-1] def append_cell(): eq_time.append([]) eq_MitoV.append([]) eq_CytoV.append([]) eq_CellV.append([]) eq_Mcl1_C.append([]) eq_Mcl1C_C.append([]) eq_Mcl1M_C.append([]) eq_Mcl1_P.append([]) eq_Mcl1C_P.append([]) eq_Mcl1M_P.append([]) result_syn_test.append([]) result_deg_test.append([]) result_time.append([]) result_MitoV.append([]) result_CytoV.append([]) result_CellV.append([]) result_Mcl1M_C.append([]) result_Mcl1C_C.append([]) result_Mcl1M_P.append([]) result_Mcl1C_P.append([]) result_kd_new.append([]) def plot_results(n): plt.figure() plt.plot(result_time[n][:],result_MitoV[n][:], label = "Mito Volume") plt.plot(result_time[n][:],result_CytoV[n][:], label = "Cyto Volume") plt.plot(result_time[n][:],result_CellV[n][:], label = "Cell Volume") plt.axvline(x = (Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10)) for p in range(Number_of_cycles-1): plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+(Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", linestyle = "--", dashes = (5,10)) plt.xlabel("time [s]") plt.ylabel("Volume [µm^3]") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Compartment volumes over time of Cell Nr. " + str(n)) plt.show() plt.figure() plt.plot(result_time[n][:],result_Mcl1M_C[n][:], label = "Mitochondrial Mcl1") plt.plot(result_time[n][:],result_Mcl1C_C[n][:], label = "Cytoplasmic Mcl1") plt.axvline(x = (Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10)) for p in range(Number_of_cycles-1): plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+(Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", linestyle = "--", dashes = (5,10)) plt.xlabel("time [s]") plt.ylabel("Concentration [nM]") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Species concentrations over time of Cell Nr. " + str(n)) plt.show() plt.figure() plt.plot(result_time[n][:],result_Mcl1M_P[n][:], label = "Mitochondrial Mcl1 particles") plt.plot(result_time[n][:],result_Mcl1C_P[n][:], label = "Cytoplasmic Mcl1 particles") plt.axvline(x = (Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10)) for p in range(Number_of_cycles-1): plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Sim_start[n]), color = "grey", linestyle = "--", dashes = (5,10)) plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)+(Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2])+Sim_start[n]), color = "black", linestyle = "--", dashes = (5,10)) plt.xlabel("time [s]") plt.ylabel("Particles") plt.legend(loc="center left", bbox_to_anchor=(1.04,0.5)) plt.title("Species particles over time of Cell Nr. " + str(n)) plt.show() # Old way of saving time traces: # From simulation, every time step is saved in the sim_g1 array and temporarily in the respective temp_*** array # The temp array always start from the original time 0 through the number of cell cycles with 100 steps per second # For each cell, a Sim_start value [s] (random point in cell cycle) is defined in the loop structure # This Sim_start value is added to every time value of the result_time array # --> if Sim_Start would be 6235, insted of 0, 0.01, 0.02 --> 6235.01, 6235.02, 6235.03 in result_time[] # For cutting down the size of output arrays, new arrays are defined # for r in range (2*np.sum(Cell_Cycle_Durations)-Sim_start[n])*100,((Number_of_cycles)*np.sum(Cell_Cycle_Durations)-Sim_start[n])*100 # means for every r between (2*64800 - 6235)*100 and (3*64800 - 6235)*100 # --> through all values between 12 336 500 and 18 816 500 is looped (6 480 000 steps) # only every thousandth element is saved to the final result array --> 6480 steps (one cell cycle, every 10 s) # New way of saving time traces: # Initially save every thousandth element in result_time array (all 10 s) # Cut the edges of the array to show full cycles in final output in the same manner # original result_time will be 6235 6245 6255 -> length is (Number_of_cycles*Cell_cycle_Durations)/10 # Cut loop should go from 1233 to 1881 # (2*64800 - 6235)/10 to (3*64800 - 6235)/10 --> def cut_output(n,a,range_end): for n in range(int(Cells_per_batch)): #print("Cut_Output n:" +str(n)) result_time_f.append([]) result_time_s_f.append([]) result_MitoV_f.append([]) result_CytoV_f.append([]) result_CellV_f.append([]) result_Mcl1M_C_f.append([]) result_Mcl1C_C_f.append([]) result_Mcl1M_P_f.append([]) result_Mcl1C_P_f.append([]) result_Mcl1T_C_f.append([]) result_Mcl1T_P_f.append([]) result_syn_test_f.append([]) result_deg_test_f.append([]) # for every element from (2x cell cycle - Sim_Start) until (simulation end - Sim_Start) for r in range(math.floor((np.sum(Cell_Cycle_Durations)+Sim_start[int(n+a*Cells_per_batch)])/10),math.floor(((Number_of_cycles-1)*np.sum(Cell_Cycle_Durations)+Sim_start[int(n+a*Cells_per_batch)])/10)): result_time_f[n].append((result_time[n][r])-(np.sum(Cell_Cycle_Durations)+Sim_start[int(n+a*Cells_per_batch)])) # sets sim_start as effective t = 0 result_time_s_f[n].append(((result_time[n][r])-(np.sum(Cell_Cycle_Durations)+Sim_start[int(n+a*Cells_per_batch)]))/3600) # similar, just as hours result_MitoV_f[n].append(result_MitoV[n][r]) result_CytoV_f[n].append(result_CytoV[n][r]) result_CellV_f[n].append(result_CellV[n][r]) result_Mcl1M_C_f[n].append(result_Mcl1M_C[n][r]) result_Mcl1C_C_f[n].append(result_Mcl1C_C[n][r]) result_Mcl1M_P_f[n].append(result_Mcl1M_P[n][r]) result_Mcl1C_P_f[n].append(result_Mcl1C_P[n][r]) result_Mcl1T_C_f[n].append((((result_Mcl1M_C[n][r]*result_MitoV[n][r])/result_CellV[n][r])+((result_Mcl1C_C[n][r]*result_CytoV[n][r])/result_CellV[n][r]))) # total mcl-1 con. (weighted average) result_Mcl1T_P_f[n].append(result_Mcl1M_P[n][r]+result_Mcl1C_P[n][r]) # total Mcl-1 particles result_syn_test_f[n].append(result_syn_test[n][r]*mcl1_syn_g1) result_deg_test_f[n].append(result_deg_test[n][r]) def plot_summary_graphs(): plt.figure() plt.figure(figsize=(3,2)) plt.plot(result_time_s_f[0][:],result_MitoV_f[0][:], label = "Mitochondria", color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[0][:],result_CytoV_f[0][:], label = "Cytoplasm", color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[0][:],result_CellV_f[0][:], label = "Cell", color = "tab:green", linewidth = 1) for p in range(Number_of_cells-1): plt.plot(result_time_s_f[p+1][:],result_MitoV_f[p+1][:], color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_CytoV_f[p+1][:], color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_CellV_f[p+1][:], color = "tab:green", linewidth = 1) #plt.plot(result_time_s_f[p][:],result_CellV_f[p][:], label = "Cytoplasm", color = "tab:green", linewidth = 1) #plt.axvline(x = (Cell_Cycle_Durations[0]/3600), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10), linewidth = 0.5) #plt.axvline(x = ((Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10), linewidth = 0.5) #for p in range(Number_of_cycles-3): # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+Cell_Cycle_Durations[0]/3600), color = "grey", linestyle = "--", dashes = (5,10), linewidth = 0.5) # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+(Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", linestyle = "--", dashes = (5,10), linewidth = 0.5) plt.xlabel("time [h]", fontsize=6) plt.ylabel("Volume [µm$^3$]", fontsize=6) plt.title("Compartment volumes", fontsize=8) plt.legend(loc="upper right",prop={'size':4}) plt.xticks(fontsize=6) plt.yticks(fontsize=6) plt.tight_layout() plt.savefig(project_name + "/Compartments_CC_" + project_name + ".png", dpi=600) plt.show plt.figure() plt.figure(figsize=(3,2)) plt.plot(result_time_s_f[0][:],result_Mcl1M_C_f[0][:], label = "Mitochondria", color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[0][:],result_Mcl1C_C_f[0][:], label = "Cytoplasm", color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[0][:],result_Mcl1T_C_f[0][:], label = "Cell", color = "tab:green", linewidth = 1) for p in range(Number_of_cells-1): plt.plot(result_time_s_f[p+1][:],result_Mcl1M_C_f[p+1][:], color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_Mcl1C_C_f[p+1][:], color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_Mcl1T_C_f[p+1][:], color = "tab:green", linewidth = 1) #plt.plot(result_time_s_f[p][:],result_Mcl1T_C_f[p][:], label = "Cytoplasm", color = "tab:green", linewidth = 1) #plt.axvline(x = (Cell_Cycle_Durations[0]/3600), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10), linewidth = 0.5) #plt.axvline(x = ((Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10), linewidth = 0.5) #for p in range(Number_of_cycles-3): # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+Cell_Cycle_Durations[0]/3600), color = "grey", linestyle = "--", dashes = (5,10), linewidth = 0.5) # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+(Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", linestyle = "--", dashes = (5,10), linewidth = 0.5) plt.xlabel("time [h]", fontsize=6) plt.ylabel("Mcl1 [nM]", fontsize=6) plt.title("Mcl1 Concentration", fontsize=8) plt.legend(loc="upper right",prop={'size':4}) plt.xticks(fontsize=6) plt.yticks(fontsize=6) plt.tight_layout() plt.savefig(project_name + "/Concentrations_CC_" + project_name + ".png", dpi=600) plt.show plt.figure() plt.figure(figsize=(3,2)) plt.plot(result_time_s_f[0][:],result_Mcl1M_P_f[0][:], label = "Mitochondria", color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[0][:],result_Mcl1C_P_f[0][:], label = "Cytoplasm", color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[0][:],result_Mcl1T_P_f[0][:], label = "Cell", color = "tab:green", linewidth = 1) for p in range(Number_of_cells-1): plt.plot(result_time_s_f[p+1][:],result_Mcl1M_P_f[p+1][:], color = "tab:blue", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_Mcl1C_P_f[p+1][:], color = "tab:orange", linewidth = 1) plt.plot(result_time_s_f[p+1][:],result_Mcl1T_P_f[p+1][:], color = "tab:green", linewidth = 1) #plt.plot(result_time_s_f[p][:],result_Mcl1T_P_f[p][:], label = "Cytoplasm", color = "tab:green", linewidth = 1) #plt.axvline(x = (Cell_Cycle_Durations[0]/3600), color = "grey", label = "G1 -> S/G2", linestyle = "--", dashes = (5,10), linewidth = 0.5) #plt.axvline(x = ((Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", label = "Mitosis", linestyle = "--", dashes = (5,10), linewidth = 0.5) #for p in range(Number_of_cycles-3): # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+Cell_Cycle_Durations[0]/3600), color = "grey", linestyle = "--", dashes = (5,10), linewidth = 0.5) # plt.axvline(x = ((p+1)*np.sum(Cell_Cycle_Durations)/3600+(Cell_Cycle_Durations[0]/3600+Cell_Cycle_Durations[1]/3600+Cell_Cycle_Durations[2]/3600)), color = "black", linestyle = "--", dashes = (5,10), linewidth = 0.5) plt.xlabel("time [h]", fontsize=6) plt.ylabel("Mcl1 [particles]", fontsize=6) plt.title("Mcl1 Particles", fontsize=8) plt.legend(loc="upper right",prop={'size':4}) plt.xticks(fontsize=6) plt.yticks(fontsize=6) plt.tight_layout() plt.savefig(project_name + "/Particles_CC_" + project_name + ".png", dpi=600) plt.show def clear_old_results(): # Create output folders try: os.makedirs(project_name + "/Arrays/") except FileExistsError: pass try: os.makedirs(project_name + "/Arrays/Sim") except FileExistsError: pass filelist = [ f for f in os.listdir(project_name + "/Arrays/Sim")] for f in filelist: os.remove(os.path.join(project_name + "/Arrays/Sim", f)) try: os.makedirs(project_name + "/Arrays/EQ") except FileExistsError: pass filelist = [ f for f in os.listdir(project_name + "/Arrays/EQ")] for f in filelist: os.remove(os.path.join(project_name + "/Arrays/EQ", f)) try: os.makedirs(project_name + "/SimInfo/") except FileExistsError: pass filelist = [ f for f in os.listdir(project_name + "/SimInfo")] for f in filelist: os.remove(os.path.join(project_name + "/SimInfo", f)) # remove existing .csv from prior runs in same folder filelist = [ f for f in os.listdir(project_name + "/Arrays/") if f.endswith(".csv") ] for f in filelist: os.remove(os.path.join(project_name + "/Arrays/", f)) filelist = [ f for f in os.listdir(project_name + "/SimInfo/") if f.endswith(".csv") ] for f in filelist: os.remove(os.path.join(project_name + "/SimInfo/", f)) def save_result_arrays(Sim_start,n,a,range_end): # create 2D eq result matrix for every cell temporarily for c in range(int(Cells_per_batch)): #print("C in save_results: " +str(c)) #print("Output number: "+str(c+a*Cells_per_batch)) results = [] results.append(eq_time[c]) results.append(eq_MitoV[c]) results.append(eq_CytoV[c]) results.append(eq_CellV[c]) results.append(eq_Mcl1_C[c]) results.append(eq_Mcl1C_C[c]) results.append(eq_Mcl1M_C[c]) results.append(eq_Mcl1_P[c]) results.append(eq_Mcl1C_P[c]) results.append(eq_Mcl1M_P[c]) results = np.array(results).T.tolist() # write results with respective headers to .csv file for each cell headers = ["Time [s]", "Mito Volume [µm^3]", "Cyto Volume [µm^3]", "Cell Volume [µm^3]", "Initial Mcl-1 [nM]", "Cytoplasmic Mcl-1 [nM]", "Mitochondrial Mcl-1 [nM]", "Initial Mcl-1[n]", "Cytoplasmic Mcl-1 [n]", "Mitochondrial Mcl-1 [n]"] with open(project_name + "/Arrays/EQ/EQ_Arrays_Cell_" + str(int(c+a*Cells_per_batch)) + "_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(headers) with open(project_name + "/Arrays/EQ/EQ_Arrays_Cell_" + str(int(c+a*Cells_per_batch)) + "_" + project_name + ".csv", "a", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerows(results) # create 2D sim result matrix for every cell temporarily for c in range(int(Cells_per_batch)): results = [] results.append(result_time_f[c]) results.append(np.divide(result_time_f[c],60)) results.append(result_time_s_f[c]) results.append(result_MitoV_f[c]) results.append(result_CytoV_f[c]) results.append(result_CellV_f[c]) results.append(result_Mcl1M_C_f[c]) results.append(result_Mcl1C_C_f[c]) results.append(result_Mcl1T_C_f[c]) results.append(result_Mcl1M_P_f[c]) results.append(result_Mcl1C_P_f[c]) results.append(result_Mcl1T_P_f[c]) results.append(result_syn_test_f[c]) results.append(result_deg_test_f[c]) results = np.array(results).T.tolist() # write results with respective headers to .csv file for each cell headers = ["Time [s]", "Time[min]", "Time [h]", "Mito Volume [µm^3]", "Cyto Volume [µm^3]", "Cell Volume [µm^3]", "Mitochondrial Mcl-1 [nM]", "Cytoplasmic Mcl-1 [nM]", "Total Mcl-1 [nM]", "Mitochondrial Mcl-1[n]", "Cytoplasmic Mcl-1 [n]", "Total Mcl-1 [n]","Mcl-1 synthesis [s^-1]","Mcl-1 degradation [s^-1]"] with open(project_name + "/Arrays/Sim/Result_Arrays_Cell_" + str(int(c+a*Cells_per_batch)) + "_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(headers) with open(project_name + "/Arrays/Sim/Result_Arrays_Cell_" + str(int(c+a*Cells_per_batch)) + "_" + project_name + ".csv", "a", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerows(results) def save_siminfo(Sim_start,n,a,range_end): # Save other important meta-data to be able to re-create the simulation with open(project_name + "/SimInfo/SimStart_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(Sim_start_temp) with open(project_name + "/SimInfo/C_Mcl1_ini_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(c_mcl1_temp) with open(project_name + "/SimInfo/V_Mito_ini_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(v_mito_temp) with open(project_name + "/SimInfo/V_Cyto_ini_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(v_cyto_temp) with open(project_name + "/SimInfo/kd_ini_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerow(kd_ini_temp) with open(project_name + "/SimInfo/Treatment_" + project_name + ".csv", "w") as f: f.write("Time of treatment;" + str(treatment_time) + "\n") f.write("CHX;" + str(chx) + "\n") f.write("CHX effect;" + str(1-chx_effect) + "\n") f.write("BTZ;" + str(btz) + "\n") f.write("BTZ effect;" + str(1-btz_effect) + "\n") # Save the current notebook display(Javascript('IPython.notebook.save_checkpoint()')) shutil.copy("01_WildType_Sim.ipynb",project_name + "/SimInfo/01_WildType_Sim.ipynb") with open(project_name + "/SimInfo/Input_Parameters_" + project_name + ".csv", "w", newline="") as f: writer = csv.writer(f, delimiter =';') writer.writerows(np.transpose(Input_all)) # Simulation parameters for Simulation in general Time_Simulation_EQ = [0, 30000, 3000000] # define simulation start, end, stepsize for EQ Cell_Cycle_Durations = [29808,33048,1944] # G1, S/G2, and M phase duration [s] Cell_Cycle = Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]+Cell_Cycle_Durations[2] Number_of_cells = 500 # number of cells to simulate Number_of_batches = 25 # in how many separate simulation batches all cells are split Cells_per_batch = Number_of_cells/Number_of_batches Number_of_cycles = 4 # number of cell cycles for every cell (2 less than defined are plotted in the end) project_name = "25_M0009_WildType_500cells" # Initial values as distributions #Mcl1_dist = np.genfromtxt("Total_Mcl1.csv", delimiter=",", skip_header=1) # read Mcl-1 intensity distribution for G1 from IF imaging # Mcl-1 turnover definitions # Degradation/Synthesis rates can be adjusted here for different cell cycle stages mcl1_deg_g1 = 1.77*2.83*10**(-4) # s^-1 from literature, multiplied by 1.77 to fit to 40 min half-life mcl1_deg_sg2 = 1.77*2.83*10**(-4) # s^-1 from literature mcl1_deg_m = 1.77*2.83*10**(-4) # s^-1 from literature mcl1_syn_g1 = 1.77*2.83*10**(-4) # s^-1 from literature mcl1_syn_sg2 = 1.77*2.83*10**(-4) # s^-1 from literature mcl1_syn_m = 1.77*2.83*10**(-4) # s^-1 from literature # Degradation might be altered over cell cycle, can be adjusted # through the following variables (they directly multiply the degradation rate) deg_g1 = 1 deg_sg2 = 1 deg_m = 0.4 # Synthesis might be altered over cell cycle, can be adjusted # through the following variables (they directly multiply the synthesis rate) syn_g1 = 1 syn_sg2 = 1 syn_m = 1 # Define treatment time, sort of treatment and it's impact treatment_time = 43200 # [s] chx = 0 # CHX treatment? 0 = no, 1 = yes btz = 0 # BTZ treatment? 0 = no, 1 = yes chx_effect = 0 # CHX effect from 0 (100 %) to 1 (0 %) btz_effect = 0 # BTZ effect from 0 (100 %) to 1 (0 %) # Define the KD value between Mito and Cyto (lower means more drive towards Mito) # adjustable for different cell cycle phases # initial kd_sg2 and kd_m are not used atm (defined by prior phases, but still need to be defined) kd_eq = 0.58#np.random.normal(loc=0.6,scale=0.2) kd_g1 = kd_eq kd_sg2 = kd_eq kd_m = kd_eq # Define the speed of Mcl1 interchanging between Cyto and Mito (until now arbitrarely defined) k_off_eq = 15 k_off = 1.23 # experimental k_off = 1.2393 --> times 2 (assuming average ratio of 0.5) --> k_on = 2.4786 # Define the rate at which KD decreases from the initial value kd_g1 in the # respective cell cycle phase (End of one phase is used as start for next phase) kd_rate_g1 = 0#3.333333333333333333333333*10**(-6) kd_rate_sg2 = 0#1.438888888888888888888888*10**(-5) kd_rate_m = 0#5.138888888888888888888888*10**(-5) # Define growth rates per volume and second for different cell cycle phases # No growth in M, value still needs to be defined growth_g1 = 9.16666666666666666666666*10**(-6) #s^-1 growth_sg2 = 1.138888888888888888888888*10**(-5) #s^-1 growth_m = 1.138888888888888888888888*10**(-5) #s^-1 # Print all input variables in csv to check quickly Input_head = ["G1_duration", "S/G2_duration", "M_duration", "Number_of_Cells", "Number_of_Cycles",\ "Mcl1_deg_G1","Mcl1_deg_S/G2","Mcl1_deg_M","Mcl1_syn_G1","Mcl1_syn_S/G2","Mcl1_syn_M",\ "Mcl1_degF_G1","Mcl1_degF_S/G2","Mcl1_degF_M","Mcl1_synF_G1","Mcl1_synF_S/G2","Mcl1_synF_M",\ "t_time","chx","btz","chx_effect","btz_effect","kd_eq","kd_g1","kd_S/G2","kd_M",\ "k_off_eq","k_off","kd_rate_G1","kd_rate_S/G2","kd_rate_M","growth_G1","growth_S/G2","growth_M"] Input_values = [Cell_Cycle_Durations[0],Cell_Cycle_Durations[1],Cell_Cycle_Durations[2],Number_of_cells,Number_of_cycles,\ mcl1_deg_g1,mcl1_deg_sg2,mcl1_deg_m,mcl1_syn_g1,mcl1_syn_sg2,mcl1_syn_m,\ deg_g1,deg_sg2,deg_m,syn_g1,syn_sg2,syn_m,\ treatment_time,chx,btz,chx_effect,btz_effect,kd_eq,kd_g1,kd_sg2,kd_m,\ k_off_eq,k_off,kd_rate_g1,kd_rate_sg2,kd_rate_m,growth_g1,growth_sg2,growth_m] Input_all = [Input_head,Input_values] def save_params(): try: os.makedirs(project_name + "/Parameters/") except FileExistsError: pass start = time.process_time() clear_old_results() c_mcl1_ini_list = [] v_mito_ini_list = [] v_cyto_ini_list = [] kd_ini_list = [] Sim_start = [] for n in range(Number_of_cells): kd_ini_list.append(np.random.normal(loc=kd_eq,scale=0.02)) c_mcl1_ini_list.append(np.random.normal(loc=110,scale=10)) # add initial mcl-1 conc. to array for saving later mito_temp = np.random.normal(loc=466,scale=54.89) v_mito_ini_list.append(mito_temp) # add initial mito vol. to array for saving later v_cyto_ini_list.append(mito_temp*7.4) # add initial cyto vol. to array for saving later Sim_start.append(10)#random.randrange(np.sum(Cell_Cycle_Durations))) # choose a random starting time point within the cell cycle duration #print(Sim_start) for a in range(Number_of_batches): # Definition of arrays that store output varibales as time, volume, concentrations # over multiple cell cycles eq_time = [] eq_MitoV = [] eq_CytoV = [] eq_CellV = [] eq_Mcl1_C = [] eq_Mcl1C_C = [] eq_Mcl1M_C = [] eq_Mcl1_P = [] eq_Mcl1C_P = [] eq_Mcl1M_P = [] result_time = [] result_MitoV = [] result_CytoV = [] result_CellV = [] result_Mcl1M_C = [] result_Mcl1C_C = [] result_Mcl1M_P = [] result_Mcl1C_P = [] result_kd_new = [] result_syn_test = [] result_deg_test = [] result_time_f = [] result_time_s_f = [] result_MitoV_f = [] result_CytoV_f = [] result_CellV_f = [] result_Mcl1M_C_f = [] result_Mcl1C_C_f = [] result_Mcl1T_C_f = [] result_Mcl1M_P_f = [] result_Mcl1C_P_f = [] result_Mcl1T_P_f = [] result_syn_test_f = [] result_deg_test_f = [] range_end = (Cells_per_batch)+a*Cells_per_batch #print("Batch Nr. (a): " +str(a)) #print("Range_end: " +str(range_end)) Sim_start_temp = Sim_start[0:int((a+1)*Cells_per_batch)] c_mcl1_temp = c_mcl1_ini_list[0:int((a+1)*Cells_per_batch)] v_mito_temp = v_mito_ini_list[0:int((a+1)*Cells_per_batch)] v_cyto_temp = v_cyto_ini_list[0:int((a+1)*Cells_per_batch)] kd_ini_temp = kd_ini_list[0:int((a+1)*Cells_per_batch)] for n in range(int(range_end-Cells_per_batch), int(range_end)): #print("Cell N: " + str(n)) # Define the KD value between Mito and Cyto (lower means more drive towards Mito) # adjustable for different cell cycle phases # initial kd_sg2 and kd_m are not used atm (defined by prior phases, but still need to be defined) kd_eq = kd_ini_list[n] kd_g1 = kd_ini_list[n] kd_sg2 = kd_ini_list[n] kd_m = kd_ini_list[n] c_mcl1_ini = c_mcl1_ini_list[n] # choose mcl-1 concentration out of the distribution v_mito_ini = v_mito_ini_list[n] # define mitochondrial volume between 420 - 512 µm^3 (error range) v_cyto_ini = v_cyto_ini_list[n] # define cytoplasmic volume as 7.4-fold (experimental ratio) #print(v_mito_ini) append_cell() # functions adds new 2D arrays to results arrays (3D) - every cell has 2D arrays within the 3D result arrays load_eq() # loads Antimony model structure of eq_simulation eq_preparation() # defines input values within the loaded Antimony model for equilibration selections = ['time'] + ['Mito', 'Cyto', 'Cell'] + toymodel_eq.getFloatingSpeciesConcentrationIds() + toymodel_eq.getFloatingSpeciesIds() # defines which parameters of the whole simulation are saved sim_eq = toymodel_eq.simulate(start=Time_Simulation_EQ[0],end=Time_Simulation_EQ[1],steps=Time_Simulation_EQ[2], selections = selections) # the model is simulated save_results_eq(int(n-a*Cells_per_batch)) # results of equilibration simulation needed for the real simulation are saved (Compartment volumes, Mcl-1 concentrations and particle numbers) t_time = np.sum(Cell_Cycle_Durations) + treatment_time + Sim_start[n] # scale the predefined treatment time, so that the treatment occurs within the second cell cycle and normalized to the respective sim_start value - so that all cells get treated at the same time time_passed = 0 # variable that counts, how much time has passed during simulation treat = 0 # switch-like variable to decide, when to start the treatment --> used in functions treat_con = 0 # switch-like variable, when to continue with treatment --> used in functions treat_dif = 0 # used later - defines the total simulation time at the start of each new phase simulation to time the treatment correctly for c in range(Number_of_cycles): # loops through cell cycles # Prepare and simulate G1 phase if (t_time >= np.sum(Cell_Cycle_Durations)) and (t_time <= (np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0])) and (time_passed >= np.sum(Cell_Cycle_Durations)) and (time_passed <= (np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0])): # if the treatment time is within the G1 phase of the second cycle AND the time passed maches the treatment range treat = 1 #treatment is activated treat_dif = np.sum(Cell_Cycle_Durations) # treatment difference to total simulation start (1 CC) is defined for later substraction if (t_time >= (2*np.sum(Cell_Cycle_Durations))) and (t_time <= (2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0])) and (time_passed >= 2*np.sum(Cell_Cycle_Durations)) and (time_passed <= (2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0])): # if the treatment time is within the G1 phase of the third cycle AND the time passed maches the treatment range treat = 1 #treatment is activated treat_dif = 2*np.sum(Cell_Cycle_Durations) # treatment difference to total simulation start (2 CC) is defined for later substraction load_g1() # loads the antimony model structure of G1 simulation g1_preparation(t_time, treat, time_passed, treat_con, treat_dif) #loads input variables for G1 simulation selections_g1 = ['time'] + ['Mito', 'Cyto', 'Cell'] + toymodel_sim_g1.getFloatingSpeciesConcentrationIds() + toymodel_sim_g1.getFloatingSpeciesIds() + toymodel_sim_g1.getGlobalParameterIds() sim_g1 = toymodel_sim_g1.simulate(start=(0),end=(Cell_Cycle_Durations[0]),steps=(Cell_Cycle_Durations[0]*100), selections = selections_g1) # model G1 phase save_results_g1(c,np.sum(Cell_Cycle_Durations),int(n-a*Cells_per_batch)) # save results/parameters from G1 time_passed += Cell_Cycle_Durations[0] # increased the time_passed by duration of G1 if treat == 1: # if treatment happened in G1 treat_con = 1 # continuing treatment is activated # Prepare and simulate S/G2 phase - similar to G1 if (t_time >= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]) and (t_time <= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]) and (time_passed >= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]) and (time_passed <= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]): treat = 1 treat_dif = np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0] if (t_time >= 2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]) and (t_time <= 2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]) and (time_passed >= 2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]) and (time_passed <= 2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]): treat = 1 treat_dif = 2*np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0] load_sg2() sg2_preparation(t_time, treat, time_passed, treat_con, treat_dif) selections_sg2 = ['time'] + ['Mito', 'Cyto', 'Cell'] + toymodel_sim_sg2.getFloatingSpeciesConcentrationIds() + toymodel_sim_sg2.getFloatingSpeciesIds() + toymodel_sim_sg2.getGlobalParameterIds() sim_sg2 = toymodel_sim_sg2.simulate(start=(0),end=(Cell_Cycle_Durations[1]),steps=(Cell_Cycle_Durations[1]*100), selections = selections_sg2) save_results_sg2(c,np.sum(Cell_Cycle_Durations),Cell_Cycle_Durations[0],int(n-a*Cells_per_batch)) time_passed += Cell_Cycle_Durations[1] if treat == 1: treat_con = 1 # Prepare and simulate M phase - similar to G1 if (t_time >= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]) and (t_time <= 2*np.sum(Cell_Cycle_Durations)) and (time_passed >= np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1]) and (time_passed <= 2*np.sum(Cell_Cycle_Durations)): treat = 1 treat_dif = np.sum(Cell_Cycle_Durations)+Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1] load_m() m_preparation(t_time, treat, time_passed, treat_con, treat_dif) selections_m = ['time'] + ['Mito', 'Cyto', 'Cell'] + toymodel_sim_m.getFloatingSpeciesConcentrationIds() + toymodel_sim_m.getFloatingSpeciesIds() + toymodel_sim_m.getGlobalParameterIds() sim_m = toymodel_sim_m.simulate(start=(0),end=(Cell_Cycle_Durations[2]),steps=(Cell_Cycle_Durations[2]*100), selections=selections_m) save_results_m(c,np.sum(Cell_Cycle_Durations),Cell_Cycle_Durations[0]+Cell_Cycle_Durations[1],int(n-a*Cells_per_batch)) time_passed += Cell_Cycle_Durations[2] if treat == 1: treat_con = 1 print("Finished Cell Nr. " + str(n) + " Elapsed Time: " + str(round((time.process_time() - start)/60, 2)) + " min") cut_output(n,a,range_end) # cut output traces --> from (2x cell cycle - Sim_Start) until (simulation end - Sim_Start) save_result_arrays(Sim_start,n,a,range_end) save_siminfo(Sim_start,n,a,range_end) print("Batch Run-time: " + str(round((time.process_time() - start)/60, 2)) + " min") print("Final Run-time: " + str(round((time.process_time() - start)/60, 2)) + " min") save_params() plt.figure() plt.figure(figsize=(3,2)) plt.plot(eq_time[0],eq_Mcl1_C[0], label = "Initial Mcl-1", color = "tab:red", linewidth = 1.5) plt.plot(eq_time[0],eq_Mcl1M_C[0], label = "Mito Mcl-1", color = "tab:blue", linewidth = 1.5) plt.plot(eq_time[0],eq_Mcl1C_C[0], label = "Cyto Mcl-1", color = "tab:orange", linewidth = 1.5) plt.xlabel("time [s]", fontsize=6) plt.ylabel("Mcl1 [nM]", fontsize=6) plt.title("Mcl1 Concentration", fontsize=8) plt.legend(loc="upper right",prop={'size':4}) plt.xticks(fontsize=6) plt.yticks(fontsize=6) plt.tight_layout() plt.savefig(project_name + "/Concentrations_EQ_" + project_name + ".png", dpi=600) plt.show