Skip to content

Commit f0aa6ea

Browse files
author
Hiroaki Imoto
authored
Release v0.6.0 (#46)
* Update README.md * Create scipy_differential_evolution * Remove visualize.jl and run_simulation() * Use isinstalled when including pyproject * Remove convert.jl * Reshape array containing simulation results (#45) * Reshape simulation array * Fix LoadError: syntax * Revert visualize.jl * Revert PyPlot * Fix dirname * Update python functions * Use println * Fix index * Don't use warn * Move isinstalled * Fix index * Create numpy_load * Remove initpop * Update README.md * Add generate_initial_population * Bump version to 0.6.0 * Remove println
1 parent a176daa commit f0aa6ea

File tree

11 files changed

+570
-431
lines changed

11 files changed

+570
-431
lines changed

‎Project.toml

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "BioMASS"
22
uuid = "324734c7-f323-4536-9335-775d9be9d101"
33
authors = ["Hiroaki Imoto <himoto@protein.osaka-u.ac.jp>"]
4-
version = "0.5.1"
4+
version = "0.6.0"
55

66
[deps]
77
CMAEvolutionStrategy = "8d3b24bd-414e-49e0-94fb-163cc3a3e411"
@@ -14,7 +14,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1414
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
1515
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
1616
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
17-
Seaborn = "d2ef9438-c967-53ab-8060-373fdd9e13eb"
1817
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1918
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
2019
SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
@@ -28,7 +27,6 @@ DelayDiffEq = "5.28"
2827
ForwardDiff = "0.10"
2928
PyCall = "1.92"
3029
PyPlot = "2.9"
31-
Seaborn = "0.4, 1"
3230
StatsBase = "0.33"
3331
SteadyStateDiffEq = "1.5"
3432
Sundials = "4.3"

‎README.md

Lines changed: 15 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,16 +21,26 @@ or through the `pkg` REPL mode by typing
2121
] add BioMASS
2222
```
2323

24+
### Python package requirements:
25+
26+
- numpy - https://numpy.org
27+
- scipy - https://scipy.org
28+
- matplotlib - https://matplotlib.org
29+
2430
## Example
2531

2632
### Model development
2733

28-
This example shows you how to build a simple Michaelis-Menten two-step enzyme catalysis model. [`pasmopy.Text2Model`](https://pasmopy.readthedocs.io/en/latest/model_development.html) allows you to build a BioMASS model from text. You simply describe biochemical reactions and the molecular mechanisms extracted from text are converted into an executable model.
34+
This example shows you how to build a simple Michaelis-Menten two-step enzyme catalysis model.
35+
36+
> E + S ⇄ ES → E + P
37+
38+
[`pasmopy.Text2Model`](https://pasmopy.readthedocs.io/en/latest/model_development.html) allows you to build a BioMASS model from text. You simply describe biochemical reactions and the molecular mechanisms extracted from text are converted into an executable model.
2939

3040
Prepare a text file describing the biochemical reactions (e.g., `michaelis_menten.txt`)
3141
```
32-
E binds S <--> ES | kf=0.003, kr=0.001 | E=100, S=50
33-
ES dissociates to E and P | kf=0.002, kr=0
42+
E + S ⇄ ES | kf=0.003, kr=0.001 | E=100, S=50
43+
ES → E + P | kf=0.002
3444
3545
@obs Substrate: u[S]
3646
@obs E_free: u[E]
@@ -72,8 +82,8 @@ using BioMASS
7282

7383
model = Model("./examples/fos_model");
7484

75-
# Estimate unknown model parameters against experimental observations.
76-
optimize(model, 1, max_generation=20000, allowable_error=0.5)
85+
# Estimate unknown model parameters from experimental observations
86+
scipy_differential_evolution(model, 1) # requires scipy package
7787

7888
# Save simulation results to figure/ in the model folder
7989
run_simulation(model, viz_type="best", show_all=true)

‎examples/fos_model/fitness.jl

Lines changed: 24 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
# Residual Sum of Squares
22
function compute_objval_rss(
3-
sim_data::Vector{Float64},
4-
exp_data::Vector{Float64})::Float64
3+
sim_data::Vector{Float64},
4+
exp_data::Vector{Float64})::Float64
55
error::Float64 = 0.0
66
for i in eachindex(exp_data)
77
@inbounds error += (sim_data[i] - exp_data[i])^2
@@ -12,9 +12,9 @@ end
1212

1313
# Cosine similarity
1414
function compute_objval_cos(
15-
sim_data::Vector{Float64},
16-
exp_data::Vector{Float64})::Float64
17-
error::Float64 = 1.0 - dot(sim_data,exp_data)/(norm(sim_data)*norm(exp_data))
15+
sim_data::Vector{Float64},
16+
exp_data::Vector{Float64})::Float64
17+
error::Float64 = 1.0 - dot(sim_data, exp_data) / (norm(sim_data) * norm(exp_data))
1818
return error
1919
end
2020

@@ -28,59 +28,59 @@ end
2828

2929

3030
function diff_sim_and_exp(
31-
sim_matrix::Matrix{Float64},
32-
exp_dict::Dict{String,Array{Float64,1}},
33-
exp_timepoint::Vector{Float64},
34-
conditions::Vector{String};
35-
sim_norm_max::Float64)::Tuple{Vector{Float64}, Vector{Float64}}
31+
sim_matrix::Matrix{Float64},
32+
exp_dict::Dict{String,Array{Float64,1}},
33+
exp_timepoint::Vector{Float64},
34+
conditions::Vector{String};
35+
sim_norm_max::Float64)::Tuple{Vector{Float64},Vector{Float64}}
3636
sim_result::Vector{Float64} = []
3737
exp_result::Vector{Float64} = []
3838

39-
for (idx,condition) in enumerate(conditions)
39+
for (idx, condition) in enumerate(conditions)
4040
if condition in keys(exp_dict)
41-
append!(sim_result,sim_matrix[Int.(exp_timepoint.+1),idx])
42-
append!(exp_result,exp_dict[condition])
41+
append!(sim_result, sim_matrix[idx, Int.(exp_timepoint .+ 1)])
42+
append!(exp_result, exp_dict[condition])
4343
end
4444
end
4545

46-
return (sim_result./sim_norm_max, exp_result)
46+
return (sim_result ./ sim_norm_max, exp_result)
4747
end
4848

4949

5050
# Define an objective function to be minimized.
51-
function objective(indiv_gene)::Float64
51+
function objective(indiv_gene)::Float64
5252
indiv::Vector{Float64} = decode_gene2val(indiv_gene)
5353

54-
(p,u0) = update_param(indiv)
54+
(p, u0) = update_param(indiv)
5555

56-
if Sim.simulate!(p,u0) isa Nothing
56+
if Sim.simulate!(p, u0) isa Nothing
5757
error::Vector{Float64} = zeros(length(observables))
58-
for (i,obs_name) in enumerate(observables)
59-
if isassigned(Exp.experiments,i)
58+
for (i, obs_name) in enumerate(observables)
59+
if isassigned(Exp.experiments, i)
6060
if length(Sim.normalization) > 0
6161
norm_max::Float64 = (
6262
Sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
6363
Sim.simulations[
6464
i,
65-
Sim.normalization[obs_name]["timepoint"],
66-
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
65+
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]],
66+
Sim.normalization[obs_name]["timepoint"]
6767
]
6868
) : maximum(
6969
Sim.simulations[
7070
i,
71+
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]],
7172
:,
72-
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
7373
]
7474
)
7575
)
7676
end
7777
error[i] = compute_objval_rss(
7878
diff_sim_and_exp(
79-
Sim.simulations[i,:,:],
79+
Sim.simulations[i, :, :],
8080
Exp.experiments[i],
8181
Exp.get_timepoint(obs_name),
8282
Sim.conditions,
83-
sim_norm_max = ifelse(
83+
sim_norm_max=ifelse(
8484
length(Sim.normalization) == 0, 1.0, norm_max
8585
)
8686
)...

‎examples/fos_model/simulation.jl

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -28,20 +28,20 @@ const t = collect(0.0:dt:5400.0) # 0, 1, 2, ..., 5400 [sec.]
2828
const conditions = ["EGF", "HRG"]
2929

3030
simulations = Array{Float64,3}(
31-
undef, length(observables), length(t), length(conditions)
31+
undef, length(observables), length(conditions), length(t)
3232
)
3333

3434

3535
function solveode(
36-
f::Function,
37-
u0::Vector{Float64},
38-
t::Vector{Float64},
39-
p::Vector{Float64})::Union{ODESolution{},Nothing}
36+
f::Function,
37+
u0::Vector{Float64},
38+
t::Vector{Float64},
39+
p::Vector{Float64})::Union{ODESolution{},Nothing}
4040
local sol::ODESolution{}, is_successful::Bool
4141
try
4242
prob = ODEProblem(f, u0, (t[1], t[end]), p)
4343
sol = solve(
44-
prob,CVODE_BDF(),
44+
prob, CVODE_BDF(),
4545
abstol=ABSTOL,
4646
reltol=RELTOL,
4747
saveat=dt,
@@ -61,11 +61,11 @@ end
6161

6262

6363
function get_steady_state(
64-
f::Function,
65-
u0::Vector{Float64},
66-
p::Vector{Float64})::Vector{Float64}
64+
f::Function,
65+
u0::Vector{Float64},
66+
p::Vector{Float64})::Vector{Float64}
6767
local sol::SteadyStateSolution{}, is_successful::Bool
68-
try
68+
try
6969
prob = ODEProblem(f, u0, (0.0, Inf), p)
7070
prob = SteadyStateProblem(prob)
7171
sol = solve(
@@ -110,29 +110,29 @@ function simulate!(p::Vector{Float64}, u0::Vector{Float64})::Union{Bool,Nothing}
110110
return false
111111
else
112112
@inbounds @simd for j in eachindex(t)
113-
simulations[observables_index("Phosphorylated_MEKc"),j,i] = (
113+
simulations[observables_index("Phosphorylated_MEKc"), i, j] = (
114114
sol.u[j][V.ppMEKc]
115115
)
116-
simulations[observables_index("Phosphorylated_ERKc"),j,i] = (
116+
simulations[observables_index("Phosphorylated_ERKc"), i, j] = (
117117
sol.u[j][V.pERKc] + sol.u[j][V.ppERKc]
118118
)
119-
simulations[observables_index("Phosphorylated_RSKw"),j,i] = (
119+
simulations[observables_index("Phosphorylated_RSKw"), i, j] = (
120120
sol.u[j][V.pRSKc] + sol.u[j][V.pRSKn] * (p[C.Vn] / p[C.Vc])
121121
)
122-
simulations[observables_index("Phosphorylated_CREBw"),j,i] = (
122+
simulations[observables_index("Phosphorylated_CREBw"), i, j] = (
123123
sol.u[j][V.pCREBn] * (p[C.Vn] / p[C.Vc])
124124
)
125-
simulations[observables_index("dusp_mRNA"),j,i] = (
125+
simulations[observables_index("dusp_mRNA"), i, j] = (
126126
sol.u[j][V.duspmRNAc]
127127
)
128-
simulations[observables_index("cfos_mRNA"),j,i] = (
128+
simulations[observables_index("cfos_mRNA"), i, j] = (
129129
sol.u[j][V.cfosmRNAc]
130130
)
131-
simulations[observables_index("cFos_Protein"),j,i] = (
131+
simulations[observables_index("cFos_Protein"), i, j] = (
132132
(sol.u[j][V.pcFOSn] + sol.u[j][V.cFOSn]) * (p[C.Vn] / p[C.Vc])
133133
+ sol.u[j][V.cFOSc] + sol.u[j][V.pcFOSc]
134134
)
135-
simulations[observables_index("Phosphorylated_cFos"),j,i] = (
135+
simulations[observables_index("Phosphorylated_cFos"), i, j] = (
136136
sol.u[j][V.pcFOSn] * (p[C.Vn] / p[C.Vc]) + sol.u[j][V.pcFOSc]
137137
)
138138
end

‎examples/nfkb_model/fitness.jl

Lines changed: 25 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ function diff_sim_and_exp(
3838

3939
for (idx, condition) in enumerate(conditions)
4040
if condition in keys(exp_dict)
41-
append!(sim_result, sim_matrix[Int.(exp_timepoint .+ 1),idx])
41+
append!(sim_result, sim_matrix[idx, Int.(exp_timepoint .+ 1)])
4242
append!(exp_result, exp_dict[condition])
4343
end
4444
end
@@ -48,7 +48,7 @@ end
4848

4949

5050
# Define an objective function to be minimized.
51-
function objective(indiv_gene)::Float64
51+
function objective(indiv_gene)::Float64
5252
indiv::Vector{Float64} = decode_gene2val(indiv_gene)
5353

5454
(p, u0) = update_param(indiv)
@@ -59,32 +59,32 @@ function objective(indiv_gene)::Float64
5959
if isassigned(Exp.experiments, i)
6060
if length(Sim.normalization) > 0
6161
norm_max::Float64 = (
62-
Sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
63-
Sim.simulations[
64-
i,
65-
Sim.normalization[obs_name]["timepoint"],
66-
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
67-
]
68-
) : maximum(
69-
Sim.simulations[
70-
i,
71-
:,
72-
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]]
73-
]
62+
Sim.normalization[obs_name]["timepoint"] !== nothing ? maximum(
63+
Sim.simulations[
64+
i,
65+
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]],
66+
Sim.normalization[obs_name]["timepoint"]
67+
]
68+
) : maximum(
69+
Sim.simulations[
70+
i,
71+
[conditions_index(c) for c in Sim.normalization[obs_name]["condition"]],
72+
:,
73+
]
74+
)
7475
)
75-
)
7676
end
7777
error[i] = compute_objval_rss(
78-
diff_sim_and_exp(
79-
Sim.simulations[i,:,:],
80-
Exp.experiments[i],
81-
Exp.get_timepoint(obs_name),
82-
Sim.conditions,
83-
sim_norm_max=ifelse(
84-
length(Sim.normalization) == 0, 1.0, norm_max
85-
)
86-
)...
87-
)
78+
diff_sim_and_exp(
79+
Sim.simulations[i, :, :],
80+
Exp.experiments[i],
81+
Exp.get_timepoint(obs_name),
82+
Sim.conditions,
83+
sim_norm_max=ifelse(
84+
length(Sim.normalization) == 0, 1.0, norm_max
85+
)
86+
)...
87+
)
8888
end
8989
end
9090
return sum(error) # < 1e12

‎examples/nfkb_model/simulation.jl

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -26,26 +26,26 @@ const sstime = 1000.0 # time to reach steady state
2626
const conditions = ["WT"]
2727

2828
simulations = Array{Float64,3}(
29-
undef, length(observables), length(t), length(conditions)
29+
undef, length(observables), length(conditions), length(t)
3030
)
3131

3232
function solvedde(
33-
f::Function,u0::Vector{Float64},history::Vector{Float64},
34-
tspan::Tuple{Float64,Float64},p::Vector{Float64},tau::Float64)
33+
f::Function, u0::Vector{Float64}, history::Vector{Float64},
34+
tspan::Tuple{Float64,Float64}, p::Vector{Float64}, tau::Float64)
3535
h(p, t) = history
3636
lags = [tau]
3737
prob = DDEProblem(f, u0, h, tspan, p; constant_lags=lags)
3838
alg = MethodOfSteps(BS3())
3939
sol = solve(
40-
prob,alg,saveat=dt,abstol=1e-8,reltol=1e-8,verbose=false
40+
prob, alg, saveat=dt, abstol=1e-8, reltol=1e-8, verbose=false
4141
)
4242
return sol
4343
end
4444

4545

4646
function get_steady_state(
47-
p::Vector{Float64},u0::Vector{Float64},
48-
sstime::Float64,tau::Float64)::Vector{Float64}
47+
p::Vector{Float64}, u0::Vector{Float64},
48+
sstime::Float64, tau::Float64)::Vector{Float64}
4949
# get steady state (t<0)
5050
p[C.term] = 1.0
5151
history::Vector{Float64} = u0
@@ -64,8 +64,8 @@ end
6464

6565

6666
function get_time_course(
67-
p::Vector{Float64},u0::Vector{Float64},
68-
sstime::Float64,tau::Float64)
67+
p::Vector{Float64}, u0::Vector{Float64},
68+
sstime::Float64, tau::Float64)
6969
p1::Vector{Float64} = copy(p)
7070
p1[C.term] = 0.0
7171
u1::Vector{Float64} = get_steady_state(p, u0, sstime, tau)
@@ -84,8 +84,8 @@ end
8484

8585

8686
function simulate!(
87-
p::Vector{Float64},
88-
u0::Vector{Float64})::Union{Bool,Nothing}
87+
p::Vector{Float64},
88+
u0::Vector{Float64})::Union{Bool,Nothing}
8989
for (i, condition) in enumerate(conditions)
9090
# if condition == "WT"
9191
# pass
@@ -95,7 +95,7 @@ function simulate!(
9595
return false
9696
else
9797
@inbounds @simd for j in eachindex(t)
98-
simulations[observables_index("Nuclear_NFkB"),j,i] = (
98+
simulations[observables_index("Nuclear_NFkB"), i, j] = (
9999
sol.u[j][V.NFKBn]
100100
)
101101
end

0 commit comments

Comments
 (0)