It’s a simulation of gene regulatory network using Gillespie’s stochastic simulation algorithm.

```
using Statistics
using Dates
using CSV
using DataFrames
function FillArrays0(Nanog,Gata6,Recep,Fgf4,idx,t)
for i in 1:1
if (i != idx)
Nanog[i,t+1] = Nanog[i,t]
Gata6[i,t+1] = Gata6[i,t]
Recep[i,t+1] = Recep[i,t]
Fgf4[i,t+1] = Fgf4[i,t]
end
end
return Nanog, Gata6, Recep, Fgf4
end
function Run0()
println("I'm in m0")
#numT = 1
#FracVec_II = Array{Float64}(undef,numT)
T = 50
#YMax = 100
Cols = 10^8 + 2*10^6
Val_kf = 96; Val_Kd2f = 25; Val_Km = 100
Vec = [1.25, 0.3, 1.25, 0.3, 2.0, 0.3, 160, 5.0, 17, 13.5, Val_kf/32, 1, Val_Kd2f, Val_Km]
p = 2
b1=p*Vec[1]
d1=Vec[2]
b2=p*Vec[3]
d2=Vec[4]
b3=p*Vec[5]
d3=Vec[6]
Kd21=p*Vec[7] #binding of Gata6 to Nanog's promoter
Kd31=p*Vec[8] #binding of Activated Recep to Nanog's promoter
Kd12=p*Vec[9] #binding of Nanog to Gata6's promoter
Kd13=p*Vec[10] #binding of Nanog to Recep's promoter
kf=p*Vec[11]
d4=Vec[12]
Kd2f=p*Vec[13] #binding of Gata6 to Fgf4's promoter
Km=p*Vec[14]
a=2 #Hill-coefficient of Nanog
b=2 #Hill-coefficient of Gata6
c=2 #Hill-coefficient of Fgf4 suppression of Nanog production
d=2 #Hill-coefficient of Gata6 repression of Fgf4 production
#name = 'Stoch_Gillespie_1b'
#filename1 = strcat('D:\Research\Excel_files\stochastic_stable_states\',name,'_stochastic_stable_states.xlsx')
#for num = 1:1:numT
#println(num)
t=1
Nanog = Array{Float64}(undef,1,Cols)
Gata6 = Array{Float64}(undef,1,Cols)
Recep = Array{Float64}(undef,1,Cols)
Fgf4 = Array{Float64}(undef,1,Cols)
Time = Array{Float64}(undef,1,Cols)
Nanog_Total = 0
Gata6_Total = 0
for i1 in 1:1
for j1 in 1:Cols
Nanog[i1,j1] = -1
Gata6[i1,j1] = -1
Recep[i1,j1] = -1
Fgf4[i1,j1] = -1
end
end
for i2 in 1:Cols
Time[1,i2] = -1
end
for i3 in 1:size(Nanog)[1]
Nanog[i3,1] = 0
Gata6[i3,1] = 0
Recep[i3,1] = 0
Fgf4[i3,1] = 0
end
Time[1,1] = 0
True = true
while True
PF = zeros(1,8)
#XSn is the average number of molecules perceived by the nth cell
XS1 = 0.5*Fgf4[1,t]
m = 0
#PF[,n) is the propensity function for the (quotient(PF[,),6)+1)th cell
Ra1 = Recep[1,t]*(XS1/Km)/(1+XS1/Km)
PF[1,1] = (2^m)*(b1)/((1+(Gata6[1,t]/Kd21)^b)*(1+(Ra1/Kd31)^c))
PF[1,2] = d1*Nanog[1,t]
PF[1,3] = (2^m)*(b2)/(1+(Nanog[1,t]/Kd12)^a)
PF[1,4] = d2*Gata6[1,t]
PF[1,5] = (2^m)*(b3)/(1+(Nanog[1,t]/Kd13)^a)
PF[1,6] = d3*Recep[1,t]
PF[1,7] = kf/(1+(Gata6[1,t]/Kd2f)^d)
PF[1,8] = d4*Fgf4[1,t]
B = rand()
C = log(B)
Tao = (-1)*C/(sum(PF))
Time[1,t+1] = Time[1,t] + Tao
diff = Time[1,t+1]-Time[1,t]
#println(Time[1,t+1])
#if (diff < small)
# small=diff
#end
k=size(Time)
if ((Time[1,t]+Tao) >= T)
#print("Time: ")
#println(Time[1,t])
#print("Tao: ")
#println(Tao)
break
end
D = rand()
if (D*sum(PF) < PF[1,1])
Nanog[1,t+1] = Nanog[1,t]+1
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (PF[1,1] <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:2]))
Nanog[1,t+1] = Nanog[1,t]-1
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (sum(PF[1,1:2]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:3]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]+1
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (sum(PF[1,1:3]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:4]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]-1
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (sum(PF[1,1:4]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:5]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]+1
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (sum(PF[1,1:5]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:6]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]-1
Fgf4[1,t+1] = Fgf4[1,t]
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
elseif (sum(PF[1,1:6]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:7]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]+1
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
else #(sum(PF[1,1:7]) <= D*sum(PF) && D*sum(PF) < sum(PF[1,1:8]))
Nanog[1,t+1] = Nanog[1,t]
Gata6[1,t+1] = Gata6[1,t]
Recep[1,t+1] = Recep[1,t]
Fgf4[1,t+1] = Fgf4[1,t]-1
Ret = FillArrays0(Nanog,Gata6,Recep,Fgf4,1,t)
Nanog = Ret[1]
Gata6 = Ret[2]
Recep = Ret[3]
Fgf4 = Ret[4]
end
t = t+1
#println(t)
end
#display(Nanog[1,1:5])
EndIdx = 0
for i4 in 1:size(Nanog)[2]
if(Nanog[1,i4]==-1)
EndIdx = i4-1
break
end
end
#print("EndIdx: ")
#println(EndIdx)
NewNanog = zeros(1,EndIdx)
NewGata6 = zeros(1,EndIdx)
NewRecep = zeros(1,EndIdx)
NewFgf4 = zeros(1,EndIdx)
NewTime = zeros(1,EndIdx)
for i5 in 1:1
for j2 in 1:EndIdx
NewNanog[i5,j2] = Nanog[i5,j2]
NewGata6[i5,j2] = Gata6[i5,j2]
NewRecep[i5,j2] = Recep[i5,j2]
NewFgf4[i5,j2] = Fgf4[i5,j2]
end
end
for j2 in 1:EndIdx
NewTime[1,j2] = Time[1,j2]
end
RetNanog_m0 = zeros(1)
RetGata6_m0 = zeros(1)
RetRecep_m0 = zeros(1)
RetFgf4_m0 = zeros(1)
for i7 in 1:1
RetNanog_m0[i7] = NewNanog[i7,end]
RetGata6_m0[i7] = NewGata6[i7,end]
RetRecep_m0[i7] = NewRecep[i7,end]
RetFgf4_m0[i7] = NewFgf4[i7,end]
end
#=
for i6 in 1:1
#println(string("Image_SS_II_",i6))
if (NewNanog[i6,EndIdx] > NewGata6[i6,EndIdx])
Nanog_Total = Nanog_Total + 1
else
Gata6_Total = Gata6_Total + 1
end
#=
plot(NewNanog[i6,:], lc=:blue, lw=2)
plot!(NewGata6[i6,:], lc=:green, lw=2)
plot!(NewRecep[i6,:], lc=:red, lw=2)
plot!(NewFgf4[i6,:], lc=:magenta ,lw=2)
plot!(grid=false, legend=false)
xlims!(0, size(NewNanog)[2])
ylims!(0, YMax)
xlabel!("Time")
ylabel!("Protein concentration")
savefig(string("Documents/Fgf4/Plots/Stoch_II/Fig_Stoch_II_NPF_Plot_",i6,".png"))
=#
#=
plot(NewNanog[i6,:], lc=:blue, lw=2)
plot!(NewGata6[i6,:], lc=:green, lw=2)
plot!(grid=false, legend=false)
xlims!(0, size(NewNanog)[2])
ylims!(0, YMax)
xlabel!("Time")
ylabel!("Protein concentration")
savefig(string("/home/Dinkar/Documents/Fgf4/Plots/Stoch_II/Fig_Stoch_II_NPF_Plot_Sans_Recep_Fgf4_",i6,".png"))
=#
end
Frac = Nanog_Total/(Nanog_Total + Gata6_Total)
FracVec_II[num] = Frac
println(Frac)
=#
return RetNanog_m0, RetGata6_m0, RetRecep_m0, RetFgf4_m0
#end
end
```

I am simulating a gene regulatory network along with cell division, so with every cell division the number of equations doubles. This is for the first cell. I have four more such scripts (for no. of cells 2, 4, 8, 16, 32). And I access all the files from the file for 32 cell stage.

Plz let me know how this code could be sped up using mutli-processing (at this point I am too scared of something going wrong with multi-threading).