{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Build the model\n", "\n", "using DifferentialEquations\n", "using DiffEqBiological\n", "using DelimitedFiles\n", "using Plots\n", "using Distributions\n", "using Statistics\n", "using XLSX\n", "using StatsPlots\n", "using KernelDensity\n", "using Random\n", "\n", "Bcell_model = @reaction_network begin\n", " \n", "#V\n", " V*hill(pMyc*(100/V),k1,Kd1,n1), ∅ --> V\n", " \n", "#S\n", " k0, S --> ∅\n", "\n", "#CTV\n", " 0, ∅ --> CTV\n", " 0, CTV --> ∅\n", "\n", "#IRF4 \n", " (V/V0)*(k10 + mm(S,k11,Kd11) + hill(pIRF4*(V0/V),k12,Kd12,n12) + hill(pPRDM1*(V0/V),k33,Kd33,n33)*hillr(pIRF8*(V0/V),k36,Kd36,n36)), ∅ --> mIRF4\n", " k13, mIRF4 --> ∅\n", " k14*mIRF4, ∅ --> pIRF4\n", " k15, pIRF4 --> ∅\n", " \n", "#IRF8 \n", " (V/V0)*(k40 + mm(S,k41,Kd41) + hillr(pIRF4*(V0/V),k42,Kd42,n42) * hill(pPAX5*(V0/V),k43,Kd43,n43)), ∅ --> mIRF8\n", " k44, mIRF8 --> ∅\n", " k45*mIRF8, ∅ --> pIRF8\n", " k46, pIRF8 --> ∅\n", " \n", "#BCL6 \n", " (V/V0)*(k16 + k17*hill(pIRF4*(V0/V),1,Kd13,n13)*hillr(pPRDM1*(V0/V),1,Kd131,n131)*hillr(pIRF4*(V0/V),1,Kd133,n133) + hill(pIRF8*(V0/V),k134,Kd134,n134)), ∅ --> mBCL6\n", " k18, mBCL6 --> ∅\n", " k19*mBCL6, ∅ --> pBCL6\n", " k20*(1+mm(BCR,v20,Kd132)), pBCL6 --> ∅ \n", "\n", "#PRDM1 \n", " (V/V0)*(k21 + k21*mm(BCR,1,Kd132) + mm(S,k35,Kd35) + hill(pIRF4*(V0/V),k32,Kd32,n32)*hillr(Division_counter*(V0/V),k202,Kd202,n202) + (hillr(pBach2*(V0/V),k34,Kd34,n34) + hillr(pPAX5*(V0/V),k22,Kd14,n14))*hillr(pBCL6*(V0/V),k31,Kd31,n31)), ∅ --> mPRDM1\n", " k23, mPRDM1 --> ∅\n", " k24*mPRDM1, ∅ --> pPRDM1\n", " k25, pPRDM1 --> ∅\n", " \n", "#PAX5\n", " (V/V0)*(k26 + hillr(pPRDM1*(V0/V),k27,Kd15,n15)), ∅ --> mPAX5\n", " k28, mPAX5 --> ∅\n", " k29*mPAX5, ∅ --> pPAX5\n", " k30, pPAX5 --> ∅\n", "\n", "#Myc\n", " (V/V0)*(k101 + hill(pBCL6*(V0/V),k102,Kd102,n102)*hillr(pPRDM1*(V0/V),k103,Kd103,n103)*hill(Division_counter*(V0/V),k202,Kd202,n202)), ∅ --> mMyc\n", " #(V/V0)*(k101 + hill(pBCL6*(V0/V),k102,Kd102,n102)*hill(Division_counter*(V0/V),k203,Kd203,n203)), ∅ --> mMyc\n", " k104, mMyc --> ∅\n", " k105*mMyc, ∅ --> pMyc\n", " k106, pMyc --> ∅ \n", "\n", "#Bach2\n", " (V/V0)*(k70 + hill(pBCL6*(V0/V),k71,Kd71,n71) + hill(pPAX5*(V0/V),k72,Kd72,n72) + hill(TCDD,k301,Kd301,n301)), ∅ --> mBach2\n", " k73, mBach2 --> ∅\n", " k74*mBach2, ∅ --> pBach2\n", " k75, pBach2 --> ∅ \n", "\n", "#AID\n", " (V/V0)*(k80 + (hill(pIRF8*(V0/V),k86,Kd86,n86) + hill(pBCL6*(V0/V),k81,Kd81,n81))*hill(pBach2*(V0/V),k82,Kd82,n82)), ∅ --> mAID\n", " k83, mAID --> ∅\n", " k84*mAID, ∅ --> pAID\n", " k85, pAID --> ∅ \n", "\n", "#IgM\n", " (V/V0)*(k90 + hillr(pPAX5*(V0/V),k91,Kd91,n91)), ∅ --> mIgM\n", " k92, mIgM --> ∅\n", " k93*mIgM, ∅ --> pIgM\n", " k94, pIgM --> ∅ \n", " \n", "#Division_counter\n", " k200, ∅ --> Division_counter\n", " k201, Division_counter --> ∅ \n", " \n", " end #=\n", " =# V0 #=\n", " =# v20 #=\n", " =# BCR #=\n", " =# TCDD #=\n", " =# k101 #=\n", " =# k102 #=\n", " =# k103 #=\n", " =# k104 #=\n", " =# k105 #=\n", " =# k106 #=\n", " =# k0 #=\n", " =# k1 #=\n", " =# k2 #=\n", " =# k10 #=\n", " =# k11 #=\n", " =# k12 #=\n", " =# k13 #=\n", " =# k14 #=\n", " =# k15 #=\n", " =# k16 #=\n", " =# k17 #=\n", " =# k18 #=\n", " =# k19 #=\n", " =# k20 #=\n", " =# k21 #=\n", " =# k22 #=\n", " =# k23 #=\n", " =# k24 #=\n", " =# k25 #=\n", " =# k26 #=\n", " =# k27 #=\n", " =# k28 #=\n", " =# k29 #=\n", " =# k30 #=\n", " =# k31 #=\n", " =# k32 #=\n", " =# k33 #=\n", " =# k34 #=\n", " =# k35 #=\n", " =# k36 #=\n", " =# k40 #=\n", " =# k41 #=\n", " =# k42 #=\n", " =# k43 #=\n", " =# k44 #=\n", " =# k45 #=\n", " =# k46 #=\n", " =# k70 #=\n", " =# k71 #=\n", " =# k72 #=\n", " =# k73 #=\n", " =# k74 #=\n", " =# k75 #=\n", " =# k80 #=\n", " =# k81 #=\n", " =# k82 #=\n", " =# k83 #=\n", " =# k84 #=\n", " =# k85 #=\n", " =# k86 #=\n", " =# k90 #=\n", " =# k91 #=\n", " =# k92 #=\n", " =# k93 #=\n", " =# k94 #=\n", " =# k134 #=\n", " =# k200 #=\n", " =# k201 #=\n", " =# k202 #=\n", " =# k203 #=\n", " =# k301 #=\n", "\n", " =# Kd1 #=\n", " =# Kd2 #=\n", " =# Kd11 #=\n", " =# Kd12 #= \n", " =# Kd13 #=\n", " =# Kd131 #=\n", " =# Kd132 #=\n", " =# Kd133 #=\n", " =# Kd134 #= \n", " =# Kd14 #=\n", " =# Kd15 #=\n", " =# Kd16 #=\n", " =# Kd31 #=\n", " =# Kd32 #=\n", " =# Kd33 #=\n", " =# Kd34 #=\n", " =# Kd35 #=\n", " =# Kd36 #=\n", " =# Kd41 #=\n", " =# Kd42 #=\n", " =# Kd43 #=\n", " =# Kd71 #=\n", " =# Kd72 #=\n", " =# Kd81 #=\n", " =# Kd82 #=\n", " =# Kd86 #=\n", " =# Kd91 #=\n", " =# Kd102 #=\n", " =# Kd103 #= \n", " =# Kd202 #=\n", " =# Kd203 #=\n", " =# Kd301 #=\n", "\n", " =# n1 #=\n", " =# n2 #=\n", " =# n12 #=\n", " =# n13 #=\n", " =# n131 #=\n", " =# n133 #=\n", " =# n134 #=\n", " =# n14 #=\n", " =# n15 #=\n", " =# n16 #=\n", " =# n31 #= \n", " =# n32 #=\n", " =# n33 #=\n", " =# n34 #=\n", " =# n36 #=\n", " =# n42 #=\n", " =# n43 #=\n", " =# n71 #=\n", " =# n72 #=\n", " =# n81 #=\n", " =# n82 #=\n", " =# n86 #=\n", " =# n91 #=\n", " =# n102 #=\n", " =# n103 #=\n", " =# n202 #=\n", " =# n203 #=\n", " =# n301 \n", "\n", "parr = Dict([ \n", " (:V, 100),\n", " (:V0, 100), \n", " (:v20, 10),\n", " (:BCR, 0),\n", " (:TCDD, 0),\n", " (:k101, 0.01*log(2)/3600),\n", " (:k102, 10*log(2)/3600), \n", " (:k103, 1), \n", " (:k104, log(2)/3600), \n", " (:k105, 10*log(2)/3600/2), \n", " (:k106, log(2)/3600/2),\n", " (:k0, log(2)/3600/240),\n", " (:k1, log(2)/3600/5),\n", " #(:k1, 0),\n", " (:k2, 1.9254e-5),\n", " (:k10, 0.01*log(2)/3600), \n", " (:k11, log(2)/3600), \n", " (:k12, 10*log(2)/3600), \n", " (:k13, log(2)/3600), \n", " (:k14, 10*log(2)/3600/2), \n", " (:k15, log(2)/3600/2), \n", " (:k16, 0.01*log(2)/3600), \n", " (:k17, 10*log(2)/3600), \n", " (:k18, log(2)/3600), \n", " (:k19, 10*log(2)/3600/2), \n", " (:k20, log(2)/3600/2), \n", " (:k21, 0.04*log(2)/3600), \n", " (:k22, 0.5), \n", " (:k23, log(2)/3600), \n", " (:k24, 10*log(2)/3600/2), \n", " (:k25, log(2)/3600/2), \n", " (:k26, 0.05*log(2)/3600), \n", " (:k27, 10*log(2)/3600), \n", " (:k28, log(2)/3600), \n", " (:k29, 10*log(2)/3600/2), \n", " (:k30, log(2)/3600/2), \n", " (:k31, 10*log(2)/3600), \n", " (:k32, 3*log(2)/3600), \n", " #(:k32, 7.7016e-7),\n", " (:k33, 10*log(2)/3600),\n", " (:k34, 0.5),\n", " (:k35, 7.7016e-7),\n", " (:k36, 1),\n", " (:k40, 2*log(2)/3600), \n", " (:k41, 12*log(2)/3600), \n", " (:k42, 1), \n", " (:k43, 5*log(2)/3600), \n", " (:k44, log(2)/3600), \n", " (:k45, 10*log(2)/3600/2), \n", " (:k46, log(2)/3600/2), \n", " (:k70, 0.01*log(2)/3600),\n", " (:k71, 10*log(2)/3600),\n", " (:k72, 10*log(2)/3600),\n", " (:k73, log(2)/3600), \n", " (:k74, 10*log(2)/3600/2), \n", " (:k75, log(2)/3600/2),\n", " (:k80, 7.7016e-7),\n", " (:k81, 7.7016e-4),\n", " (:k82, 1),\n", " (:k83, 7.7016e-5), \n", " (:k84, 1.9254e-3), \n", " (:k85, 1.9254e-4),\n", " (:k86, 7.7016e-4),\n", " (:k90, 3.8508e-7), \n", " (:k91, 3.8508e-4), \n", " (:k92, 3.8508e-5), \n", " (:k93, 3.8508e-4), \n", " (:k94, 3.8508e-5), \n", " (:k134, 1*log(2)/3600),\n", " (:k200, log(2)/3600),\n", " (:k201, log(2)/3600/10000),\n", " (:k202, 1),\n", " (:k203, 1),\n", " (:k301, 10*log(2)/3600),\n", " \n", " (:Kd1, 25),\n", " (:Kd2, 25),\n", " (:Kd11, 200),\n", " (:Kd12, 25), \n", " (:Kd13, 25), \n", " (:Kd131, 25), \n", " (:Kd132, 25), \n", " (:Kd133, 300),\n", " (:Kd134, 200),\n", " (:Kd14, 25), \n", " (:Kd15, 25), \n", " (:Kd16, 15),\n", " (:Kd31, 25),\n", " (:Kd32, 25),\n", " (:Kd33, 50),\n", " (:Kd34, 25),\n", " (:Kd35, 200),\n", " (:Kd36, 25),\n", " (:Kd41, 200),\n", " (:Kd42, 300),\n", " (:Kd43, 25),\n", " (:Kd71, 25), \n", " (:Kd72, 25),\n", " (:Kd81, 25), \n", " (:Kd82, 25),\n", " (:Kd86, 200),\n", " (:Kd91, 25),\n", " (:Kd102, 25), \n", " (:Kd103, 25),\n", " (:Kd202, 120), #in PRDM1, by 8th division(9th generation), Division_counter should drop to 20000/256 = 78\n", " (:Kd203, 110), #in Myc, by 8th division(9th generation), Division_counter should drop to 20000/256 = 78\n", " (:Kd301, 1),\n", " \n", " (:n1, 5),\n", " (:n2, 5),\n", " (:n12, 5), \n", " (:n13, 5), \n", " (:n131, 5), \n", " (:n133, 5),\n", " (:n134, 5),\n", " (:n14, 5), \n", " (:n15, 5), \n", " (:n16, 5),\n", " (:n31, 5),\n", " (:n32, 5),\n", " (:n33, 5),\n", " (:n34, 5),\n", " (:n36, 5),\n", " (:n42, 5),\n", " (:n43, 5),\n", " (:n71, 5), \n", " (:n72, 5),\n", " (:n81, 5), \n", " (:n82, 5),\n", " (:n86, 5),\n", " (:n91, 6),\n", " (:n102, 5),\n", " (:n103, 5),\n", " (:n202, 6),\n", " (:n203, 6),\n", " (:n301, 1),\n", " ])\n", "\n", "\n", "ics = Dict([\n", " (:V, 100.0),\n", " (:S, 0.),\n", " (:CTV, 10000.0),\n", " (:mIRF4, 0), \n", " (:pIRF4, 1.0),\n", " (:mBCL6, 0), \n", " (:pBCL6, 1.0),\n", " (:mMyc, 0), \n", " (:pMyc, 1.0),\n", " (:mPRDM1, 0), \n", " (:pPRDM1, 1.0),\n", " (:mPAX5, 10.0), \n", " (:pPAX5, 100.0),\n", " (:mBach2, 10.0), \n", " (:pBach2, 100.0),\n", " (:mAID, 0), \n", " (:pAID, 1.0),\n", " (:mIgM, 0), \n", " (:pIgM, 1.0),\n", " (:Division_counter, 20000.0) \n", " ])\n", "\n", "\n", "\n", "t0 = 0. #decimal point is nessissary\n", "u0 = Array([get(ics, k, 0.0) for k=keys(speciesmap(Bcell_model))])\n", "\n", "p = Array([parr[k] for k=keys(paramsmap(Bcell_model))])\n", "\n", "#Stochastic solver\n", "prob = DiscreteProblem(Bcell_model, u0, (t0,3600*24*3), p)\n", "jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", "sol = solve(jump_prob,SSAStepper());\n", "\n", "\n", "#Determinstic solver\n", "#prob = ODEProblem(Bcell_model,u0,(t0,3600*24*3),p)\n", "#sol = solve(prob,AutoTsit5(Rosenbrock23()),reltol=1e-8,abstol=1e-8)\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#display ODEs and parameters in the model\n", "\n", "using Latexify\n", "\n", "latexify(Bcell_model, cdot=true) #display ODEs\n", "latexify(Bcell_model.params, p) #display parameter names and values\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#1 Run to the steady state\n", "\n", "Random.seed!(Int64(floor(time())))\n", "num_variables = length(sol.u[1])\n", "start_cell_num = 20\n", "LPS = 50\n", "TCDD_dose = 0.1 # nM\n", "Kd300 = 1 # nM\n", "\n", "\n", "#randome starting/initial volume distribution\n", "m_vol = 100\n", "v_vol = 12.5\n", "mu_vol = log((m_vol^2)/sqrt(v_vol+m_vol^2))\n", "sigma_vol = sqrt(log(v_vol/(m_vol^2)+1))\n", "\n", "\n", "#Initial running to basal steady state\n", "\n", "ics[:S] = 0 #No stimulation\n", "\n", "t0 = 0. #decimal point is nessissary\n", "u0 = Array([get(ics, k, 0.0) for k=keys(speciesmap(Bcell_model))])\n", "\n", "\n", "initial_time_duration = 3600*24*1\n", "\n", "u0_array = Array{Float64,2}(undef, num_variables, start_cell_num)\n", "\n", "for cell_num in 1:start_cell_num\n", " \n", " \n", " ics[:V] = round(rand(LogNormal(mu_vol,sigma_vol)))\n", " ics[:CTV] = 10000*ics[:V]/m_vol #CTV dye loading proportional to cell volumn\n", " u0 = Array([get(ics, k, 0.0) for k=keys(speciesmap(Bcell_model))])\n", " \n", " #Stochastic solver\n", " prob = DiscreteProblem(Bcell_model, u0, (t0,initial_time_duration), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol = solve(jump_prob,SSAStepper())\n", " \n", " for n in 1:num_variables\n", " global u0_array[n,cell_num] = getindex(sol.u,length(sol.u))[n] #only contains the value at end \n", " end\n", "\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#2 Run the first division\n", "\n", "time_block=3600*8\n", "end_time = 3600*24*4\n", "\n", "parr[:TCDD] = TCDD_dose\n", "p = Array([parr[k] for k=keys(paramsmap(Bcell_model))])\n", "\n", "#randome division volume distribution\n", "m_vol = 200\n", "v_vol = 50\n", "mu_vol = log((m_vol^2)/sqrt(v_vol+m_vol^2))\n", "sigma_vol = sqrt(log(v_vol/(m_vol^2)+1))\n", "\n", "cell_gen_1=fill(Array{Float64,2}(undef,0,0), (1,start_cell_num)) \n", "\n", "#p = Array([parr[k] for k=keys(paramsmap(Bcell_model))])\n", "\n", "for gen1_cell_num in 1:start_cell_num\n", " \n", " u0_array[2,gen1_cell_num] = LPS*20*(Kd300/(TCDD_dose+Kd300)) # Assign initial value of S = LPS\n", " u0 = u0_array[:,gen1_cell_num]\n", " \n", " #Stochastic solver\n", " prob = DiscreteProblem(Bcell_model, u0, (0.,time_block), p) \n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol = solve(jump_prob,SSAStepper())\n", " \n", " cell_gen_1[1,gen1_cell_num] = hcat(sol.t)\n", " for n in 1:num_variables \n", " global cell_gen_1[1,gen1_cell_num] = hcat(cell_gen_1[1,gen1_cell_num], getindex.(sol.u, n))\n", " end\n", "\n", "\n", "\n", " #N=1, 1st generation\n", "\n", " i=1 #i is time point\n", "\n", " vol_division = round(rand(LogNormal(mu_vol,sigma_vol)))\n", "\n", " while i <= length(cell_gen_1[1,gen1_cell_num][:,1]) #scan through time points \n", " if cell_gen_1[1,gen1_cell_num][i,2] < vol_division #if the cell division biomarker is less than the threshold for the current time point, no division\n", " if i == length(cell_gen_1[1,gen1_cell_num][:,1]) #if the current time point is the last time point, still no division, then rerun the model for one new time block using the last time point values as the initial values\n", " # t0 = cell_gen_1[1,gen1_cell_num][end,1]\n", " u0 = cell_gen_1[1,gen1_cell_num][end,2:end] # cell fate marker isn't added yet, so use 2:end not 2:end-1\n", " prob = DiscreteProblem(Bcell_model, u0, (0.,time_block), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol1 = solve(jump_prob,SSAStepper())\n", "\n", " if cell_gen_1[1,gen1_cell_num][end,1] + sol1.t[end] < end_time #if after running new time block, the total time is still less than end time, then append the new cell data to previous one, and go back to the while loop\n", "\n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(cell_gen_1[1,gen1_cell_num][end,1] .+ sol1.t[2:end])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:end], n))\n", " end\n", "\n", " global cell_gen_1[1,gen1_cell_num] = vcat(cell_gen_1[1,gen1_cell_num],cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", "\n", " else #if after running new time block, the total time is greater than end time, then append the new cell data up to the end time to previous one, and go back to the while loop\n", "\n", " for k= 1:length(sol1.t)\n", "\n", " if cell_gen_1[1,gen1_cell_num][end,1] + sol1.t[k] >= end_time\n", "\n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(cell_gen_1[1,gen1_cell_num][end,1] .+ sol1.t[2:k])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:k], n))\n", " end\n", "\n", " global cell_gen_1[1,gen1_cell_num] = vcat(cell_gen_1[1,gen1_cell_num],cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", "\n", " break #end For loop\n", " end\n", "\n", " end\n", " end \n", " \n", " if (cell_gen_1[1,gen1_cell_num][end,2] < vol_division) & (cell_gen_1[1,gen1_cell_num][end,1] >= end_time) #Those cells whose last time point < end time or volumn is >= vol_division will go through the while loop again \n", "\n", " #To add cell fate, fate = 1.0, no division, no death\n", " fate=fill(1.0, (length(cell_gen_1[1,gen1_cell_num][:,1]),1)) \n", " global cell_gen_1[1,gen1_cell_num]=hcat(cell_gen_1[1,gen1_cell_num], fate)\n", " break #Break from the while loop\n", " end \n", " \n", " end\n", " else\n", "\n", " global cell_gen_1[1,gen1_cell_num] = cell_gen_1[1,gen1_cell_num][1:i,:] # only save the data up to the point of division\n", "\n", " #To add cell fate, fate = 2.0, division, no death\n", " fate=fill(1.0, (length(cell_gen_1[1,gen1_cell_num][1:(end-1),1]),1)) \n", " fate=vcat(fate, 2.0)\n", " global cell_gen_1[1,gen1_cell_num]=hcat(cell_gen_1[1,gen1_cell_num], fate)\n", "\n", "\n", " break\n", " end\n", "\n", " i=i+1\n", "\n", " end\n", "end\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#3 Looping through multiple generation from 2nd generation\n", "\n", "tot_gen = 10 #Total number of generations\n", "\n", "\n", "for gen_number = 2:tot_gen\n", " \n", " #Scan through cells in last generation, and record the index for those cells that are marked as dividing\n", " global dividing_cell_index=Vector{Float64}\n", " global dividing_cell_index=zeros(0)\n", " \n", " if gen_number == 2\n", "\n", " m=length(cell_gen_1[1,:]) #the number of cells in the previous generation\n", " \n", " global mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_1[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 3\n", " \n", " m=length(cell_gen_2[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_2[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 4\n", " \n", " m=length(cell_gen_3[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_3[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 5\n", " \n", " m=length(cell_gen_4[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_4[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 6\n", " \n", " m=length(cell_gen_5[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_5[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 7\n", " \n", " m=length(cell_gen_6[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_6[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 8\n", " \n", " m=length(cell_gen_7[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_7[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 9\n", " \n", " m=length(cell_gen_8[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_8[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " \n", " elseif gen_number == 10\n", " \n", " m=length(cell_gen_9[1,:]) #the number of cells in the previous generation\n", " \n", " mother_cell_array = fill(Array{Float64,1}(undef,0), (1,m)) \n", " \n", " #Go through each cell in the last generation, take the last time point values from each cell. \n", " for cell_index = 1:m \n", " global mother_cell_array[1,cell_index] = cell_gen_9[1,cell_index][end,:]\n", " #Store the index of those cells that are dividing\n", " if (mother_cell_array[1,cell_index][end] == 2.0) & (mother_cell_array[1,cell_index][1] < end_time) #fate=2.0 for division\n", " global dividing_cell_index = vcat(dividing_cell_index, cell_index)\n", " end\n", " end\n", " end \n", " \n", "\n", "\n", "n=length(dividing_cell_index) #number of dividing cells in last generation\n", "\n", "\n", "daughter_cell_array = fill(Array{Float64,2}(undef,0,0), (2,2*n)) #daughter cell array\n", "\n", "\n", "#The random volumes spliting parameters. \n", "m2 = 0.5\n", "v2 = 0.0002\n", "mu2 = log((m2^2)/sqrt(v2+m2^2))\n", "sigma2 = sqrt(log(v2/(m2^2)+1))\n", "\n", "\n", "for j=1:n #cell index for previous generation\n", " \n", " IJulia.clear_output(false)\n", " display(gen_number)\n", " display(j)\n", " \n", " #To create the initial values for 2 daughter cells of a parent cell in the previous generation \n", " u_2a = Vector{Int64}\n", " u_2b = Vector{Int64}\n", "\n", " u_2a=zeros(0)\n", " u_2b=zeros(0)\n", "\n", " p2 = rand(LogNormal(mu2,sigma2))\n", " \n", " dividing_cell_index_number = Int(dividing_cell_index[j])\n", " for h = 2:(length(mother_cell_array[1,dividing_cell_index_number])-1) #the first column is time, the last one is cell fate, variables are in between \n", " if h == 3 #the third column is S, it is not splited, inherit the same values from parients\n", " u_2a = vcat(u_2a, mother_cell_array[1,dividing_cell_index_number][h] )\n", " u_2b = vcat(u_2b, mother_cell_array[1,dividing_cell_index_number][h] )\n", " else\n", " rand_variable_a = rand(Binomial(Int(mother_cell_array[1,dividing_cell_index_number][h]),p2)) \n", " u_2a = vcat(u_2a, rand_variable_a)\n", " u_2b = vcat(u_2b, mother_cell_array[1,dividing_cell_index_number][h] .- rand_variable_a)\n", " end\n", " end\n", " \n", " t0 = mother_cell_array[1,dividing_cell_index_number][1]\n", " \n", " #Running daughter cell_a\n", " prob = DiscreteProblem(Bcell_model, u_2a, (t0,t0+time_block), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol_a = solve(jump_prob,SSAStepper()) \n", " \n", " j_a = 2*j-1 #the cell index of the daughter_cell_a in current generation \n", " #global daughter_cell_array[1,j_a] = hcat(sol_a.t, getindex.(sol_a.u, 1),getindex.(sol_a.u, 2),getindex.(sol_a.u, 3))\n", " \n", " #Extract simulation data from sol_a\n", " daughter_cell_array[1,j_a] = hcat(sol_a.t)\n", " for n in 1:num_variables\n", " global daughter_cell_array[1,j_a] = hcat(daughter_cell_array[1,j_a], getindex.(sol_a.u, n))\n", " end\n", " \n", " i=1 #time index\n", " vol_division = round(rand(LogNormal(mu_vol,sigma_vol)))\n", " \n", " while i <= length(daughter_cell_array[1,j_a][:,1]) \n", " if daughter_cell_array[1,j_a][i,2] < vol_division #division time/threshold\n", " if i == length(daughter_cell_array[1,j_a][:,1])\n", " #t0 = daughter_cell_array[1,j_a][end,1]\n", " u0 = daughter_cell_array[1,j_a][end,2:end] # cell fate marker isn't added yet, so use 2:end not 2:end-1\n", " prob = DiscreteProblem(Bcell_model, u0, (0.,time_block), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol1 = solve(jump_prob,SSAStepper())\n", " \n", " if daughter_cell_array[1,j_a][end,1] + sol1.t[end] < end_time #if after running new time block, the total time is still less than end time, then append the new cell data to previous one.\n", " # global daughter_cell_array[1,j_a] = vcat(daughter_cell_array[1,j_a],hcat(daughter_cell_array[1,j_a][end,1] .+ sol1.t[2:end], getindex.(sol1.u[2:end], 1),getindex.(sol1.u[2:end], 2), getindex.(sol1.u[2:end], 3)))\n", " \n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(daughter_cell_array[1,j_a][end,1] .+ sol1.t[2:end])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:end], n))\n", " end\n", " global daughter_cell_array[1,j_a] = vcat(daughter_cell_array[1,j_a], cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", " \n", " else \n", "\n", " for k= 1:length(sol1.t)\n", "\n", " if daughter_cell_array[1,j_a][end,1] + sol1.t[k] >= end_time\n", " # global daughter_cell_array[1,j_a] = vcat(daughter_cell_array[1,j_a], hcat(daughter_cell_array[1,j_a][end,1] .+ sol1.t[2:k], getindex.(sol1.u[2:k], 1),getindex.(sol1.u[2:k], 2), getindex.(sol1.u[2:k], 3)))\n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(daughter_cell_array[1,j_a][end,1] .+ sol1.t[2:k])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:k], n))\n", " end\n", " \n", " global daughter_cell_array[1,j_a] = vcat(daughter_cell_array[1,j_a],cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", "\n", " break #end For loop\n", " end\n", "\n", " end\n", " end \n", " \n", " if (daughter_cell_array[1,j_a][end,2] < vol_division) & (daughter_cell_array[1,j_a][end,1] >= end_time) #Those cells whose last time point < end time or volume is >= vol_division will go through the while loop again \n", " #To add cell fate, fate = 1.0, no division, no death\n", " fate=fill(1.0, (length(daughter_cell_array[1,j_a][:,1]),1)) \n", " global daughter_cell_array[1,j_a]=hcat(daughter_cell_array[1,j_a], fate)\n", " break #Break from the while loop\n", " end\n", " \n", "\n", " end\n", " \n", " else #if volumn >= vol_division\n", " \n", " global daughter_cell_array[1,j_a] = daughter_cell_array[1,j_a][1:i,:] # only save the data up to the point of division\n", " \n", " #To add cell fate, fate = 2.0, division, no death\n", " fate=fill(1.0, (length(daughter_cell_array[1,j_a][1:(end-1),1]),1)) \n", " fate=vcat(fate, 2.0)\n", " global daughter_cell_array[1,j_a]=hcat(daughter_cell_array[1,j_a], fate)\n", "\n", " break #break from the while loop\n", " end\n", "\n", " i=i+1\n", "\n", " end\n", "\n", " #To add the index of mother cell. It has to be in an array format defined by \"fill()\"\n", " global daughter_cell_array[2,j_a] = fill(dividing_cell_index_number, (1,1)) \n", " \n", " \n", " #Running daughter cell_b\n", " prob = DiscreteProblem(Bcell_model, u_2b, (t0,t0+time_block), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol_b = solve(jump_prob,SSAStepper())\n", " j_b = 2*j #the cell index of the daughter_cell_b in current generation \n", " #global daughter_cell_array[1,j_b] = hcat(sol_b.t, getindex.(sol_b.u, 1),getindex.(sol_b.u, 2), getindex.(sol_b.u, 3))\n", " #Extract simulation data from sol_b\n", " daughter_cell_array[1,j_b] = hcat(sol_b.t)\n", " for n in 1:num_variables\n", " global daughter_cell_array[1,j_b] = hcat(daughter_cell_array[1,j_b], getindex.(sol_b.u, n))\n", " end \n", " \n", " i=1 #time index\n", " vol_division = round(rand(LogNormal(mu_vol,sigma_vol)))\n", " \n", " while i <= length(daughter_cell_array[1,j_b][:,1]) \n", " if daughter_cell_array[1,j_b][i,2] < vol_division #division time/threshold\n", " if i == length(daughter_cell_array[1,j_b][:,1])\n", " # t0 = daughter_cell_array[1,j_b][end,1]\n", " u0 = daughter_cell_array[1,j_b][end,2:end] # cell fate marker isn't added yet, so use 2:end not 2:end-1\n", " prob = DiscreteProblem(Bcell_model, u0, (0.,time_block), p)\n", " jump_prob = JumpProblem(prob,Direct(),Bcell_model)\n", " sol1 = solve(jump_prob,SSAStepper())\n", " \n", " if daughter_cell_array[1,j_b][end,1] + sol1.t[end] < end_time #if after running new time block, the total time is still less than end time, then append the new cell data to previous one.\n", " #global daughter_cell_array[1,j_b] = vcat(daughter_cell_array[1,j_b],hcat(daughter_cell_array[1,j_b][end,1] .+ sol1.t[2:end], getindex.(sol1.u[2:end], 1),getindex.(sol1.u[2:end], 2), getindex.(sol1.u[2:end], 3)))\n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(daughter_cell_array[1,j_b][end,1] .+ sol1.t[2:end])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:end], n))\n", " end\n", " global daughter_cell_array[1,j_b] = vcat(daughter_cell_array[1,j_b], cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", " \n", " else \n", "\n", " for k= 1:length(sol1.t)\n", "\n", " if daughter_cell_array[1,j_b][end,1] + sol1.t[k] >= end_time\n", " #global daughter_cell_array[1,j_b] = vcat(daughter_cell_array[1,j_b], hcat(daughter_cell_array[1,j_b][end,1] .+ sol1.t[2:k], getindex.(sol1.u[2:k], 1),getindex.(sol1.u[2:k], 2), getindex.(sol1.u[2:k], 3)))\n", " \n", " #Extract simulation data from sol1 with updated time\n", " cell_gen_sol1=fill(Array{Float64,2}(undef,0,0), (1,1)) \n", " cell_gen_sol1[1,1] = hcat(daughter_cell_array[1,j_b][end,1] .+ sol1.t[2:k])\n", " for n in 1:num_variables\n", " global cell_gen_sol1[1,1] = hcat(cell_gen_sol1[1,1], getindex.(sol1.u[2:k], n))\n", " end\n", " \n", " global daughter_cell_array[1,j_b] = vcat(daughter_cell_array[1,j_b],cell_gen_sol1[1,1]) #Append sol1 data to previous run\n", " \n", " break #end For loop\n", " end\n", "\n", " end\n", " end \n", " \n", " if (daughter_cell_array[1,j_b][end,2] < vol_division) & (daughter_cell_array[1,j_b][end,1] >= end_time) #Those cells whose last time point < end time or volumn is >= vol_division will go through the while loop again \n", " #To add cell fate, fate = 1.0, no division, no death\n", " fate=fill(1.0, (length(daughter_cell_array[1,j_b][:,1]),1)) \n", " global daughter_cell_array[1,j_b]=hcat(daughter_cell_array[1,j_b], fate)\n", " break #Break from the while loop\n", " end\n", " \n", "\n", " end\n", " else\n", " \n", " global daughter_cell_array[1,j_b] = daughter_cell_array[1,j_b][1:i,:] # only save the data up to the point of division\n", " \n", " #To add cell fate, fate = 2.0, division, no death\n", " fate=fill(1.0, (length(daughter_cell_array[1,j_b][1:(end-1),1]),1)) \n", " fate=vcat(fate, 2.0)\n", " global daughter_cell_array[1,j_b]=hcat(daughter_cell_array[1,j_b], fate)\n", "\n", " break\n", " end\n", "\n", " i=i+1\n", "\n", " end\n", "\n", " global daughter_cell_array[2,j_b] = fill(dividing_cell_index_number, (1,1)) #To add the index of mother cell. It has to be in an array format defined by \"fill()\"\n", " \n", " \n", "\n", "end #end j for loop \n", " \n", " \n", " if gen_number == 2\n", " global cell_gen_2 = daughter_cell_array\n", " elseif gen_number == 3 \n", " global cell_gen_3 = daughter_cell_array\n", " elseif gen_number == 4 \n", " global cell_gen_4 = daughter_cell_array \n", " elseif gen_number == 5 \n", " global cell_gen_5 = daughter_cell_array\n", " elseif gen_number == 6 \n", " global cell_gen_6 = daughter_cell_array\n", " elseif gen_number == 7 \n", " global cell_gen_7 = daughter_cell_array\n", " elseif gen_number == 8\n", " global cell_gen_8 = daughter_cell_array\n", " elseif gen_number == 9 \n", " global cell_gen_9 = daughter_cell_array\n", " elseif gen_number == 10 \n", " global cell_gen_10 = daughter_cell_array\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#time-course plot\n", "\n", "qq =2 #variable index for plotting\n", "\n", "plot( xlabel = \"Hour\", ylabel=species(Bcell_model)[qq-1], leg=false)\n", "\n", "for vt in 1:length(cell_gen_1[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " Plots.display(plot!(cell_gen_1[1,vt][:,1]/3600, cell_gen_1[1,vt][:,qq]))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_2[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_2[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_1[1,parent_num][end,1],cell_gen_2[1,vt][:,1])/3600, vcat(cell_gen_1[1,parent_num][end,qq],cell_gen_2[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_3[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_3[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_2[1,parent_num][end,1],cell_gen_3[1,vt][:,1])/3600, vcat(cell_gen_2[1,parent_num][end,qq],cell_gen_3[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_4[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_4[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_3[1,parent_num][end,1],cell_gen_4[1,vt][:,1])/3600, vcat(cell_gen_3[1,parent_num][end,qq],cell_gen_4[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_5[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_5[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_4[1,parent_num][end,1],cell_gen_5[1,vt][:,1])/3600, vcat(cell_gen_4[1,parent_num][end,qq],cell_gen_5[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_6[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_6[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_5[1,parent_num][end,1],cell_gen_6[1,vt][:,1])/3600, vcat(cell_gen_5[1,parent_num][end,qq],cell_gen_6[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_7[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_7[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_6[1,parent_num][end,1],cell_gen_7[1,vt][:,1])/3600, vcat(cell_gen_6[1,parent_num][end,qq],cell_gen_7[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_8[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_8[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_7[1,parent_num][end,1],cell_gen_8[1,vt][:,1])/3600, vcat(cell_gen_7[1,parent_num][end,qq],cell_gen_8[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_9[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_9[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_8[1,parent_num][end,1],cell_gen_9[1,vt][:,1])/3600, vcat(cell_gen_8[1,parent_num][end,qq],cell_gen_9[1,vt][:,qq])))\n", " \n", "end \n", "\n", "for vt in 1:length(cell_gen_10[1,:])\n", " \n", " IJulia.clear_output(true)\n", " \n", " parent_num = Int(cell_gen_10[2,vt][1])\n", " \n", " Plots.display(plot!(vcat(cell_gen_9[1,parent_num][end,1],cell_gen_10[1,vt][:,1])/3600, vcat(cell_gen_9[1,parent_num][end,qq],cell_gen_10[1,vt][:,qq])))\n", " \n", "end \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Violin plot by cell generation\n", "\n", "using StatsPlots\n", "\n", "time_point_plotted =96*3600\n", "qq = 21\n", "\n", "\n", "plot(legend=false)\n", "\n", "for n in 1:10 #tot_gen\n", " \n", " global Violin_plot_data = fill(time_point_plotted, (num_variables+1))\n", " \n", " cell_generation_data = @eval ($(Symbol(\"cell_gen_$n\")))\n", " \n", " for i in 1:length(cell_generation_data[1,:])\n", " \n", " j = findfirst(j1 -> j1>=time_point_plotted, cell_generation_data[1,i][:,1]) #find the first time index where the time point >= time_point_plotted\n", " \n", " if !isnothing(j) #if j exists and \n", " if j>1 #j > 1 is used to exclude cells generated after time_point_plotted\n", " global Violin_plot_data = hcat(Violin_plot_data, cell_generation_data[1,i][j,2:end])\n", " end\n", " end \n", " end\n", " \n", " Violin_plot_data1 = Violin_plot_data[:,2:end] #remove the first randomly filled column \n", " \n", " if !isempty(Violin_plot_data1)\n", "\n", " \n", " x = [\"Generation \"*string(n)]\n", " y = log2.(1 .+ Violin_plot_data1[qq,:])\n", " \n", " #y = Violin_plot_data1[qq,:]\n", " \n", " IJulia.clear_output(true) \n", " forStructPlot = violin!(x, y, side=:both, ylabel=\"Log2 (Copy Number+1) \"*string(species(Bcell_model)[qq]), marker=(0.2,:blue,stroke(0)))\n", " dotplot!(x, y, marker=(:black,stroke(0)))\n", " \n", " Plots.display(forStructPlot)\n", " end\n", " \n", "end\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Calculate cell numbers at indicated sampling time points\n", "\n", "time_select_interval = 6*3600\n", "number_sampling_time_point = cld(end_time,time_select_interval)+1 #cld(x,y) is for rounding up (x/y)\n", "Ig_threhold = 20.\n", "\n", "B = fill(0, (number_sampling_time_point))\n", "A_total = fill(0, (number_sampling_time_point))\n", "\n", " \n", " for n in 1:10\n", " cell_generation_data = @eval ($(Symbol(\"cell_gen_$n\"))) \n", "\n", " for num in 1:length(cell_generation_data[1,:])\n", " \n", " for t = 1:number_sampling_time_point\n", " if cell_generation_data[1,num][1,1] <= (t-1)*time_select_interval # The first time point of a cell has to be either before or equal to the selected sampling time point. If the first time point of the cell is behind the selected sampling time point, then the cell deos not belong to this sampling time point.\n", " global j = findfirst(j1 -> j1>=(t-1)*time_select_interval, cell_generation_data[1,num][:,1]) #find the first time index where the time point >= time_point_plotted\n", "\n", " if !isnothing(j) #if j exists\n", " \n", " global A_total[t] = A_total[t] + 1\n", " if cell_generation_data[1,num][j,22] >= Ig_threhold\n", " global B[t] = B[t]+1\n", " end\n", " \n", " end\n", " end\n", " end\n", "\n", " end\n", "\n", " end\n", "\n", "\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#Export TXT\n", "\n", "using DelimitedFiles\n", "\n", "save_interval = 100\n", "\n", "cell_separator = fill(\"*\", (1,24))\n", "\n", "\n", "for n in 1:10\n", " cell_generation_data = @eval ($(Symbol(\"cell_gen_$n\"))) \n", " \n", " \n", " A = fill(0, (1,24))\n", " num_of_cells = length(cell_generation_data[1,:])\n", " \n", " for num in 1:num_of_cells\n", " \n", " if n==1\n", " global mother_num = fill(0, (1,24)) \n", " else\n", " global mother_num = fill(getindex(Int64.(cell_generation_data[2,num])),(1,24))\n", " end\n", " \n", " A_temp = Int64.(floor.(cell_generation_data[1,num])) \n", " j = 0\n", " \n", " \n", " for i in 2:(length(A_temp[:,1])-1) #keep the first and the last time point\n", " if i % save_interval != 0\n", " j = vcat(j,i) \n", " end\n", " end\n", " \n", " A_temp = A_temp[setdiff(1:end, j), :]\n", " \n", " A_temp = vcat(cell_separator, vcat(mother_num, A_temp))\n", " \n", " global A = vcat(A, A_temp) \n", " \n", " IJulia.clear_output(false)\n", " display(n)\n", " display(num)\n", " \n", " end\n", " \n", " global A = A[2:end,:] #remove the first row of 0s\n", " writedlm(\"cell_gen_\"*string(n)*\".txt\", A) \n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#retrive data from text\n", "\n", "cell_gen_10=[]\n", "\n", "#cell_separator = fill(999, (1,24))\n", "using DelimitedFiles\n", "\n", "for n in 1:10\n", " \n", " global cell_txt = readdlm(\"cell_gen_\"*string(n)*\".txt\")\n", " \n", " j = 0 #mother_num \n", " \n", " for i in 1:length(cell_txt[:,1])\n", " if cell_txt[i,end] == \"*\"\n", " global j = vcat(j, i+1) \n", " end\n", " end\n", " \n", " global j = vcat(j[2:end], length(cell_txt[:,1])+2)\n", " \n", " global cell_gen_txt_temp = fill(Array{Float64,2}(undef,0,0), (2,length(j)-1)) \n", " \n", " for num in 1:(length(j)-1)\n", "\n", " global cell_gen_txt_temp[2,num] = fill(cell_txt[j[num],1], (1,1))\n", " global cell_gen_txt_temp[1,num] = cell_txt[(j[num]+1):(j[num+1]-2),:] \n", " end\n", " \n", " @eval ($(Symbol(\"cell_gen_$n\"))) = cell_gen_txt_temp\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.2.0", "language": "julia", "name": "julia-1.2" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 2 }