From 33fb40fb0fe9bf19ed0a554c5a6e720a15ec6f0d Mon Sep 17 00:00:00 2001 From: Hiroaki Imoto Date: Wed, 29 Jun 2022 00:33:36 +0900 Subject: [PATCH] From text to model (#136) * Migrate Text2Model * Create test_text2model * Create __init__.py * Fix tests * Update assert * Update test for cfos_model * Migrate text2model docs * Add construction to index * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- biomass/__init__.py | 1 + biomass/construction/__init__.py | 1 + biomass/construction/julia_template.py | 543 ++++++ biomass/construction/reaction_rules.py | 1596 +++++++++++++++++ biomass/construction/template/README.md | 15 + biomass/construction/template/__init__.py | 5 + .../template/name2idx/__init__.py | 2 + .../template/name2idx/parameters.py | 19 + .../construction/template/name2idx/species.py | 19 + biomass/construction/template/observable.py | 78 + biomass/construction/template/ode.py | 34 + biomass/construction/template/problem.py | 103 ++ .../construction/template/reaction_network.py | 39 + biomass/construction/template/search_param.py | 78 + biomass/construction/template/viz.py | 79 + biomass/construction/text2model.py | 904 ++++++++++ .../thermodynamic_restrictions.py | 205 +++ docs/api/construction.rst | 10 + docs/api/index.rst | 1 + docs/api/reaction_rules.rst | 7 + docs/api/template.rst | 20 + docs/api/text2model.rst | 6 + docs/api/thermodynamic_restrictions.rst | 6 + tests/test_text2model/C.py | 235 +++ tests/test_text2model/__init__.py | 0 .../simulations_original_BN.npy | Bin 0 -> 20456 bytes tests/test_text2model/test_error_message.py | 57 + .../test_fos_model_from_txt.py | 107 ++ .../test_model_construction.py | 184 ++ .../test_thermodynamic_restriction.py | 54 + .../text_files/Kholodenko1999.txt | 44 + tests/test_text2model/text_files/abc.txt | 4 + .../text_files/duplicate_binding.txt | 2 + .../test_text2model/text_files/fos_model.txt | 80 + .../text_files/michaelis_menten.txt | 10 + tests/test_text2model/text_files/typo_1.txt | 1 + tests/test_text2model/text_files/typo_2.txt | 1 + 37 files changed, 4550 insertions(+) create mode 100644 biomass/construction/__init__.py create mode 100644 biomass/construction/julia_template.py create mode 100644 biomass/construction/reaction_rules.py create mode 100644 biomass/construction/template/README.md create mode 100644 biomass/construction/template/__init__.py create mode 100644 biomass/construction/template/name2idx/__init__.py create mode 100644 biomass/construction/template/name2idx/parameters.py create mode 100644 biomass/construction/template/name2idx/species.py create mode 100644 biomass/construction/template/observable.py create mode 100644 biomass/construction/template/ode.py create mode 100644 biomass/construction/template/problem.py create mode 100644 biomass/construction/template/reaction_network.py create mode 100644 biomass/construction/template/search_param.py create mode 100644 biomass/construction/template/viz.py create mode 100644 biomass/construction/text2model.py create mode 100644 biomass/construction/thermodynamic_restrictions.py create mode 100644 docs/api/construction.rst create mode 100644 docs/api/reaction_rules.rst create mode 100644 docs/api/template.rst create mode 100644 docs/api/text2model.rst create mode 100644 docs/api/thermodynamic_restrictions.rst create mode 100755 tests/test_text2model/C.py create mode 100644 tests/test_text2model/__init__.py create mode 100644 tests/test_text2model/simulations_original_BN.npy create mode 100644 tests/test_text2model/test_error_message.py create mode 100644 tests/test_text2model/test_fos_model_from_txt.py create mode 100644 tests/test_text2model/test_model_construction.py create mode 100644 tests/test_text2model/test_thermodynamic_restriction.py create mode 100644 tests/test_text2model/text_files/Kholodenko1999.txt create mode 100644 tests/test_text2model/text_files/abc.txt create mode 100644 tests/test_text2model/text_files/duplicate_binding.txt create mode 100644 tests/test_text2model/text_files/fos_model.txt create mode 100644 tests/test_text2model/text_files/michaelis_menten.txt create mode 100644 tests/test_text2model/text_files/typo_1.txt create mode 100644 tests/test_text2model/text_files/typo_2.txt diff --git a/biomass/__init__.py b/biomass/__init__.py index 14f3862f..d6b7cbc4 100644 --- a/biomass/__init__.py +++ b/biomass/__init__.py @@ -1,5 +1,6 @@ """BioMASS - A Python Framework for Modeling and Analysis of Signaling Systems""" +from .construction import Text2Model from .core import * from .result import OptimizationResults from .version import __version__ diff --git a/biomass/construction/__init__.py b/biomass/construction/__init__.py new file mode 100644 index 00000000..0c10b9c2 --- /dev/null +++ b/biomass/construction/__init__.py @@ -0,0 +1 @@ +from .text2model import Text2Model diff --git a/biomass/construction/julia_template.py b/biomass/construction/julia_template.py new file mode 100644 index 00000000..cf4ff592 --- /dev/null +++ b/biomass/construction/julia_template.py @@ -0,0 +1,543 @@ +""" +Model template for conversion of biomass models into +BioMASS.jl (https://github.com/biomass-dev/BioMASS.jl) format. + +BioMASS.jl >= 0.7.2 +""" + +from typing import Final + +PARAMETERS: Final[ + str +] = """\ +module C + +const NAMES = [] + +for (idx,name) in enumerate(NAMES) + eval(Meta.parse("const $name = $idx")) +end + +const NUM = length(NAMES) + +end # module +""" + + +SPECIES: Final[ + str +] = """\ +module V + +const NAMES = [] + +for (idx, name) in enumerate(NAMES) + eval(Meta.parse("const $name = $idx")) +end + +const NUM = length(NAMES) + +end # module +""" + + +ODE: Final[ + str +] = """\ +function diffeq!(du, u, p, t) + v = Dict{Int64,Float64}() + + for i in 1:V.NUM + @inbounds du[i] = 0.0 + end + +end + + +function param_values()::Vector{Float64} + p::Vector{Float64} = ones(C.NUM) + + return p +end + + +function initial_values()::Vector{Float64} + u0::Vector{Float64} = zeros(V.NUM) + + return u0 +end +""" + + +OBSERVABLE: Final[ + str +] = """\ +const observables = [] + +function observables_index(observable_name::String)::Int + if !(observable_name in observables) + error("$observable_name is not defined in observables.") + end + return findfirst(isequal(observable_name),observables) +end +""" + + +SIMULATION: Final[ + str +] = """\ +module Sim +include("./name2idx/parameters.jl") +include("./name2idx/species.jl") +include("./ode.jl") +include("./observable.jl") + +using .C +using .V + +using Sundials +using SteadyStateDiffEq + +# Options for ODE solver +const ABSTOL = 1e-8 +const RELTOL = 1e-8 + +normalization = Dict{String,Dict{}}() + +const dt = 1.0 +const t = collect(0.0:dt:100.0) + +const conditions = [] + +simulations = Array{Float64,3}( + undef, length(observables), length(conditions), length(t) +) + + +function solveode( + f::Function, + u0::Vector{Float64}, + t::Vector{Float64}, + p::Vector{Float64})::Union{ODESolution{},Nothing} + local sol::ODESolution{}, is_successful::Bool + try + prob = ODEProblem(f, u0, (t[1], t[end]), p) + sol = solve( + prob,CVODE_BDF(), + abstol=ABSTOL, + reltol=RELTOL, + saveat=dt, + dtmin=eps(), + verbose=false + ) + is_successful = ifelse(sol.retcode === :Success, true, false) + catch + is_successful = false + finally + if !is_successful + GC.gc() + end + end + return is_successful ? sol : nothing +end + + +function get_steady_state( + f::Function, + u0::Vector{Float64}, + p::Vector{Float64})::Vector{Float64} + local sol::SteadyStateSolution{}, is_successful::Bool + try + prob = ODEProblem(f, u0, (0.0, Inf), p) + prob = SteadyStateProblem(prob) + sol = solve( + prob, + DynamicSS( + CVODE_BDF(); + abstol=ABSTOL, + reltol=RELTOL + ), + dt=dt, + dtmin=eps(), + verbose=false + ) + is_successful = ifelse(sol.retcode === :Success, true, false) + catch + is_successful = false + finally + if !is_successful + GC.gc() + end + end + return is_successful ? sol.u : [] +end + + +function simulate!(p::Vector{Float64}, u0::Vector{Float64})::Union{Bool,Nothing} + # unperturbed steady state + + # add ligand + for (i, condition) in enumerate(conditions) + + sol = solveode(diffeq!, u0, t, p) + if sol === nothing + return false + else + @inbounds @simd for j in eachindex(t) + # line_num + 4 + end + end + end +end +end # module +""" + + +EXPERIMENTAL_DATA: Final[ + str +] = """\ +module Exp +include("./observable.jl") + +experiments = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables)) +error_bars = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables)) + + +function get_timepoint(obs_name::String)::Vector{Float64} + return [] +end +end # module +""" + + +SEARCH_PARAM: Final[ + str +] = """\ +# Specify model parameters and/or initial values to optimize +function get_search_index()::Tuple{Array{Int64,1},Array{Int64,1}} + # parameters + search_idx_params::Vector{Int} = [] + + # initial values + search_idx_initials::Vector{Int} = [] + + return search_idx_params, search_idx_initials +end + + +function get_search_region()::Matrix{Float64} + p::Vector{Float64} = param_values() + u0::Vector{Float64} = initial_values() + + search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index() + search_param::Vector{Float64} = initialize_search_param(search_idx, p, u0) + + search_rgn::Matrix{Float64} = zeros(2, length(p) + length(u0)) + + # Default: 0.1 ~ 10x + for (i, j) in enumerate(search_idx[1]) + search_rgn[1,j] = search_param[i] * 0.1 # lower bound + search_rgn[2,j] = search_param[i] * 10.0 # upper bound + end + + # Default: 0.5 ~ 2x + for (i, j) in enumerate(search_idx[2]) + search_rgn[1,j + length(p)] = search_param[i + length(search_idx[1])] * 0.5 # lower bound + search_rgn[2,j + length(p)] = search_param[i + length(search_idx[1])] * 2.0 # upper bound + end + + # search_rgn[:, C.param_name] = [lower_bound, upper_bound] + # search_rgn[:, V.var_name+length(p)] = [lower_bound, upper_bound] + + search_rgn = convert_scale(search_rgn, search_idx) + + return search_rgn +end + + +function update_param(indiv::Vector{Float64})::Tuple{Array{Float64,1},Array{Float64,1}} + p::Vector{Float64} = param_values() + u0::Vector{Float64} = initial_values() + + search_idx::Tuple{Array{Int64,1},Array{Int64,1}} = get_search_index() + + for (i, j) in enumerate(search_idx[1]) + @inbounds p[j] = indiv[i] + end + for (i, j) in enumerate(search_idx[2]) + @inbounds u0[j] = indiv[i + length(search_idx[1])] + end + + # parameter constraints + + return p, u0 +end + + +function decode_gene2val(indiv_gene)::Vector{Float64} + search_rgn::Matrix{Float64} = get_search_region() + indiv::Vector{Float64} = zeros(length(indiv_gene)) + + for i in eachindex(indiv_gene) + indiv[i] = 10^( + indiv_gene[i] * ( + search_rgn[2,i] - search_rgn[1,i] + ) + search_rgn[1,i] + ) + end + + return round.(indiv, sigdigits=7) +end + + +function encode_val2gene(indiv::Vector{Float64}) + search_rgn::Matrix{Float64} = get_search_region() + indiv_gene::Vector{Float64} = zeros(length(indiv)) + + for i in eachindex(indiv) + indiv_gene[i] = ( + log10(indiv[i]) - search_rgn[1,i] + ) / ( + search_rgn[2,i] - search_rgn[1,i] + ) + end + + return indiv_gene +end + + +function encode_bestIndivVal2randGene( + gene_idx::Int64, + best_indiv::Vector{Float64}, + p0_bounds::Vector{Float64})::Float64 + search_rgn::Matrix{Float64} = get_search_region() + rand_gene::Float64 = ( + log10( + best_indiv[gene_idx] * 10^( + rand() * log10(p0_bounds[2] / p0_bounds[1]) + log10(p0_bounds[1]) + ) + ) - search_rgn[1,gene_idx] + ) / ( + search_rgn[2,gene_idx] - search_rgn[1,gene_idx] + ) + return rand_gene +end + + +function initialize_search_param( + search_idx::Tuple{Array{Int64,1},Array{Int64,1}}, + p::Vector{Float64}, + u0::Vector{Float64})::Vector{Float64} + duplicate::Vector{String} = [] + if length(search_idx[1]) != length(unique(search_idx[1])) + for idx in findall( + [count(x -> x == i, search_idx[1]) for i in unique(search_idx[1])] .!= 1 + ) + push!(duplicate, C.NAMES[search_idx[1][idx]]) + end + error( + "Duplicate parameters (C.): $duplicate" + ) + elseif length(search_idx[2]) != length(unique(search_idx[2])) + for idx in findall( + [count(x -> x == i, search_idx[2]) for i in unique(search_idx[2])] .!= 1 + ) + push!(duplicate, V.NAMES[search_idx[2][idx]]) + end + error( + "Duplicate initial conditions (V.): $duplicate" + ) + end + search_param = zeros( + length(search_idx[1]) + length(search_idx[2]) + ) + for (i, j) in enumerate(search_idx[1]) + @inbounds search_param[i] = p[j] + end + for (i, j) in enumerate(search_idx[2]) + @inbounds search_param[i + length(search_idx[1])] = u0[j] + end + + if any(x -> x == 0.0, search_param) + msg::String = "search_param must not contain zero." + for idx in search_idx[1] + if p[idx] == 0.0 + error( + @sprintf( + "`C.%s` in search_idx_params: ", C.NAMES[idx] + ) * msg + ) + end + end + for idx in search_idx[2] + if u0[idx] == 0.0 + error( + @sprintf( + "`V.%s` in search_idx_initials: ", V.NAMES[idx] + ) * msg + ) + end + end + end + + return search_param +end + + +function convert_scale( + search_rgn::Matrix{Float64}, + search_idx::Tuple{Array{Int64,1},Array{Int64,1}})::Matrix{Float64} + for i = 1:size(search_rgn, 2) + if minimum(search_rgn[:,i]) < 0.0 + msg = "search_rgn[lower_bound,upper_bound] must be positive.\n" + if i <= C.NUM + error(@sprintf("`C.%s` ", C.NAMES[i]) * msg) + else + error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg) + end + elseif minimum(search_rgn[:,i]) == 0.0 && maximum(search_rgn[:,i]) != 0.0 + msg = "lower_bound must be larger than 0.\n" + if i <= C.NUM + error(@sprintf("`C.%s` ", C.NAMES[i]) * msg) + else + error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg) + end + elseif search_rgn[2,i] - search_rgn[1,i] < 0.0 + msg = "lower_bound must be smaller than upper_bound.\n" + if i <= C.NUM + error(@sprintf("`C.%s` ", C.NAMES[i]) * msg) + else + error(@sprintf("`V.%s` ", V.NAMES[i - C.NUM]) * msg) + end + end + end + + nonzero_idx::Vector{Int} = [] + for i = 1:size(search_rgn, 2) + if search_rgn[:,i] != [0.0,0.0] + push!(nonzero_idx, i) + end + end + difference::Vector{Int} = collect( + symdiff( + Set(nonzero_idx), + Set(append!(search_idx[1], C.NUM .+ search_idx[2])) + ) + ) + if length(difference) > 0 + for idx in difference + if idx <= C.NUM + println(@sprintf("`C.%s`", C.NAMES[Int(idx)])) + else + println(@sprintf("`V.%s`", V.NAMES[Int(idx) - C.NUM])) + end + end + error( + "Set these search_params in both search_idx and search_rgn." + ) + end + + search_rgn = search_rgn[:,nonzero_idx] + + return log10.(search_rgn) +end +""" + +PROBLEM: Final[ + str +] = """\ +# Residual Sum of Squares +function compute_objval_rss( + sim_data::Vector{Float64}, + exp_data::Vector{Float64})::Float64 + error::Float64 = 0.0 + for i in eachindex(exp_data) + @inbounds error += (sim_data[i] - exp_data[i])^2 + end + return error +end + + +# Cosine similarity +function compute_objval_cos( + sim_data::Vector{Float64}, + exp_data::Vector{Float64})::Float64 + error::Float64 = 1.0 - dot(sim_data, exp_data) / (norm(sim_data) * norm(exp_data)) + return error +end + + +function conditions_index(condition_name::String)::Int + if !(condition_name in Sim.conditions) + error("$condition_name is not defined in Sim.conditions") + end + return findfirst(isequal(condition_name), Sim.conditions) +end + + +function diff_sim_and_exp( + sim_matrix::Matrix{Float64}, + exp_dict::Dict{String,Array{Float64,1}}, + exp_timepoint::Vector{Float64}, + conditions::Vector{String}; + sim_norm_max::Float64)::Tuple{Vector{Float64},Vector{Float64}} + sim_result::Vector{Float64} = [] + exp_result::Vector{Float64} = [] + + for (idx, condition) in enumerate(conditions) + if condition in keys(exp_dict) + append!(sim_result, sim_matrix[idx, Int.(exp_timepoint .+ 1)]) + append!(exp_result, exp_dict[condition]) + end + end + + return (sim_result ./ sim_norm_max, exp_result) +end + + +# Define an objective function to be minimized. +function objective(indiv_gene)::Float64 + indiv::Vector{Float64} = decode_gene2val(indiv_gene) + + (p, u0) = update_param(indiv) + + if Sim.simulate!(p, u0) isa Nothing + error::Vector{Float64} = zeros(length(observables)) + for (i, obs_name) in enumerate(observables) + if isassigned(Exp.experiments, i) + if length(Sim.normalization) > 0 + norm_max::Float64 = ( + Sim.normalization[obs_name]["timepoint"] !== nothing ? maximum( + Sim.simulations[ + i, + [conditions_index(c) for c in Sim.normalization[obs_name]["condition"]], + Sim.normalization[obs_name]["timepoint"] + ] + ) : maximum( + Sim.simulations[ + i, + [conditions_index(c) for c in Sim.normalization[obs_name]["condition"]], + :, + ] + ) + ) + end + error[i] = compute_objval_rss( + diff_sim_and_exp( + Sim.simulations[i, :, :], + Exp.experiments[i], + Exp.get_timepoint(obs_name), + Sim.conditions, + sim_norm_max=ifelse( + length(Sim.normalization) == 0, 1.0, norm_max + ) + )... + ) + end + end + return sum(error) # < 1e12 + else + return 1e12 + end +end +""" diff --git a/biomass/construction/reaction_rules.py b/biomass/construction/reaction_rules.py new file mode 100644 index 00000000..abd03b37 --- /dev/null +++ b/biomass/construction/reaction_rules.py @@ -0,0 +1,1596 @@ +import re +import sys +from dataclasses import dataclass, field +from difflib import SequenceMatcher +from typing import Dict, List, NamedTuple, NoReturn, Optional + +import numpy as np + +from .thermodynamic_restrictions import ComplexFormation, DuplicateError, ThermodynamicRestrictions + + +class UnregisteredRule(NamedTuple): + expected: Optional[str] + original: Optional[str] + + +class DetectionError(Exception): + pass + + +class ArrowError(ValueError): + pass + + +PREPOSITIONS: List[str] = [ + "to", + "for", + "from", + "up", + "down", + "in", + "on", + "at", + "off", + "into", + "around", + "among", + "between", + "of", + "over", + "above", + "below", + "under", + "through", + "across", + "along", + "near", + "by", + "beside", +] + + +@dataclass +class ReactionRules(ThermodynamicRestrictions): + """Create an executable biochemical model from text. + + .. list-table:: Available reaction rules + :widths: 25 50 25 + :header-rows: 1 + + * - Rule + - Example sentence + - Parameters (optional) + * - :func:`~pasmopy.construction.reaction_rules.dimerize` + - *A* dimerizes <--> *AA* + - .. math:: kf, kr + * - :func:`~pasmopy.construction.reaction_rules.bind` + - *A* binds *B* <--> *AB* + - .. math:: kf, kr + * - :func:`~pasmopy.construction.reaction_rules.dissociate` + - *AB* dissociates to *A* and *B* + - .. math:: kf, kr + * - :func:`~pasmopy.construction.reaction_rules.is_phosphorylated` + - *uA* is phosphorylated <--> *pA* + - .. math:: kf, kr + * - :func:`~pasmopy.construction.reaction_rules.is_dephosphorylated` + - *pA* is dephosphorylated --> *uA* + - .. math:: V, K + * - :func:`~pasmopy.construction.reaction_rules.phosphorylate` + - *B* phosphorylates *uA* --> *pA* + - .. math:: V, K + * - :func:`~pasmopy.construction.reaction_rules.dephosphorylate` + - *B* dephosphorylates *pA* --> *uA* + - .. math:: V, K + * - :func:`~pasmopy.construction.reaction_rules.transcribe` + - *B* transcribes *a* + - .. math:: V, K, n, (KF, nF) + * - :func:`~pasmopy.construction.reaction_rules.synthesize` + - *B* synthesizes *A* + - .. math:: kf + * - :func:`~pasmopy.construction.reaction_rules.is_synthesized` + - *A* is synthesized + - .. math:: kf + * - :func:`~pasmopy.construction.reaction_rules.degrade` + - *B* degrades *A* + - .. math:: kf + * - :func:`~pasmopy.construction.reaction_rules.is_degraded` + - *A* is degraded + - .. math:: kf + * - :func:`~pasmopy.construction.reaction_rules.translocate` + - *Acyt* translocates from cytoplasm to nucleus (Vcyt, Vnuc) <--> *Anuc* + - .. math:: kf, kr, (V_{pre}, V_{post}) + * - :func:`~pasmopy.construction.reaction_rules.user_defined` + - @rxn Reactant --> Product: *define rate equation here* + - *user-defined* + + From v0.2.2, you can specify directionality in binding-dissociation reaction via different arrows: + + .. code-block:: python + + E + S ⇄ ES | kf=0.003, kr=0.001 | E=100, S=50 # bi-directional + ES → E + P | kf=0.002 # unidirectional + + Attributes + ---------- + input_txt : str + Model description file (*.txt), e.g., + `Kholodenko1999.txt `_. + similarity_threshold : float + If all match_scores are below this value, expected_word will not be returned. + parameters : list of strings + ``x`` : model parameters. + species : list of strings + ``y`` : model species. + reactions : list of strings + ``v`` : flux vector. + differential_equations : list of strings + ``dydt`` : right-hand side of the differential equation. + obs_desc : list of List[str] + Description of observables. + param_info : list of strings + Information about parameter values. + init_info : list of strings + Information about initial values. + param_constraints : list of strings + Information about parameter constraints. + param_excluded : list of strings + List of parameters excluded from search params because of parameter constraints. + fixed_species : list of strings + List of species which should be held fixed (never consumed) during simulation. + sim_tspan : list of strings ['t0', 'tf'] + Interval of integration. + sim_conditions : list of List[str] + Simulation conditions with stimulation. + sim_unperturbed : str + Untreated conditions to get steady state. + rule_words : dict + Words to identify reaction rules. + nothing : List[str] + Available symbol for degradation/creation to/from nothing. + fwd_arrows : List[str] + Available arrows for unidirectional reactions. + double_arrows : List[str] + Available arrows for bi-directional reactions. + + """ + + input_txt: str + similarity_threshold: float + + parameters: List[str] = field( + default_factory=list, + init=False, + ) + species: List[str] = field( + default_factory=list, + init=False, + ) + reactions: List[str] = field( + default_factory=list, + init=False, + ) + differential_equations: List[str] = field( + default_factory=list, + init=False, + ) + obs_desc: List[List[str]] = field( + default_factory=list, + init=False, + ) + param_info: List[str] = field( + default_factory=list, + init=False, + ) + init_info: List[str] = field( + default_factory=list, + init=False, + ) + param_constraints: List[str] = field( + default_factory=list, + init=False, + ) + param_excluded: List[str] = field( + default_factory=list, + init=False, + ) + fixed_species: List[str] = field( + default_factory=list, + init=False, + ) + # Information about simulation + sim_tspan: List[str] = field( + default_factory=list, + init=False, + ) + sim_conditions: List[List[str]] = field( + default_factory=list, + init=False, + ) + sim_unperturbed: str = field( + default_factory=str, + init=False, + ) + # Words to identify reaction rules + rule_words: Dict[str, List[str]] = field( + default_factory=lambda: dict( + _bind_and_dissociate=[ + " +", + ], + dimerize=[ + " dimerizes", + " homodimerizes", + " forms a dimer", + " forms dimers", + ], + bind=[ + " binds", + " forms complexes with", + ], + dissociate=[ + " is dissociated into", + " dissociates to", + ], + is_phosphorylated=[ + " is phosphorylated", + ], + is_dephosphorylated=[ + " is dephosphorylated", + ], + phosphorylate=[ + " phosphorylates", + ], + dephosphorylate=[ + " dephosphorylates", + ], + transcribe=[ + " transcribe", + " transcribes", + ], + synthesize=[ + " synthesizes", + " promotes synthesis of", + " increases", + " is translated into", + ], + is_synthesized=[ + " is synthesized", + ], + degrade=[ + " degrades", + " promotes degradation of", + " decreases", + ], + is_degraded=[ + " is degraded", + ], + translocate=[ + " translocates", + " is translocated", + ], + ), + init=False, + ) + nothing: List[str] = field( + default_factory=lambda: ["∅", "0"], + init=False, + ) + fwd_arrows: List[str] = field( + default_factory=lambda: [ + " → ", + " ↣ ", + " ↦ ", + " ⇾ ", + " ⟶ ", + " ⟼ ", + " ⥟ ", + " ⥟ ", + " ⇀ ", + " ⇁ ", + " ⇒ ", + " ⟾ ", + ] + + [" {}> ".format("-" * i) for i in range(3)], + init=False, + ) + double_arrows: List[str] = field( + default_factory=lambda: [ + " ↔ ", + " ⟷ ", + " ⇄ ", + " ⇆ ", + " ⇌ ", + " ⇋ ", + " ⇔ ", + " ⟺ ", + ] + + [" <{}> ".format("-" * i) for i in range(3)], + init=False, + ) + + def __post_init__(self) -> None: + if not 0.0 < self.similarity_threshold < 1.0: + raise ValueError("similarity_threshold must lie within (0, 1).") + + @staticmethod + def _isfloat(string: str) -> bool: + """ + Checking if a string can be converted to float. + """ + try: + float(string) + return True + except ValueError: + return False + + @staticmethod + def _remove_prefix(text: str, prefix: str) -> str: + """ + Remove prefix from a text. + """ + if text.startswith(prefix): + return text[len(prefix) :] + assert False + + def _available_arrows(self) -> List[str]: + """ + Return all available arrow types. + """ + return self.fwd_arrows + self.double_arrows + + def _set_params(self, line_num: Optional[int], func_name: Optional[str], *args: str) -> None: + """ + Set model parameters. + """ + for p_name in args: + p_name_used = p_name + ( + f"{line_num:d}" + if not p_name[-1].isdecimal() + and isinstance(line_num, int) + and func_name != "user_defined" + else "" + ) + if p_name_used not in self.parameters: + self.parameters.append(p_name_used) + + def _set_species(self, *args: str) -> None: + """ + Set model species. + """ + for s_name in args: + if s_name not in self.species: + self.species.append(s_name) + + def _raise_detection_error(self, line_num: int, line: str) -> NoReturn: + """ + Raise `DetectionError` when a keyword is invalid. + """ + unregistered_rule = self._get_partial_similarity(line) + raise DetectionError( + f"Unregistered words in line{line_num:d}: {line}" + + ( + f"\nMaybe: '{unregistered_rule.expected.lstrip()}'." + if unregistered_rule.expected is not None + else "" + ) + ) + + def _process_pval_section(self, func_name: str, line_num: int, line: str, *args: str) -> None: + + param_values = line.split("|")[1].strip().split(",") + if all("=" in pval for pval in param_values): + for pval in param_values: + base_param = pval.split("=")[0].strip(" ") + if base_param.startswith("const "): + # Parameter names with 'const' will be added to param_excluded. + base_param = base_param.split("const ")[-1] + fixed = True + else: + fixed = False + if base_param in args: + if self._isfloat(pval.split("=")[1].strip(" ")): + self.param_info.append( + "x[C." + + base_param + + (f"{line_num:d}]" if func_name != "user_defined" else "]") + + " = " + + pval.split("=")[1].strip(" ") + ) + # If a parameter value is initialized to 0.0 or fixed, + # then add it to param_excluded. + if float(pval.split("=")[1].strip(" ")) == 0.0 or fixed: + self.param_excluded.append( + base_param + + (f"{line_num:d}" if func_name != "user_defined" else "") + ) + else: + raise ValueError( + f"line{line_num:d}: Parameter value must be int or float." + ) + else: + raise ValueError( + f"line{line_num:d}: '{pval.split('=')[0].strip(' ')}'\n" + f"Available parameters are: {', '.join(args)}." + ) + elif param_values[0].strip(" ").isdecimal(): + # Parameter constraints + for param_name in args: + # base_pname = self._get_base_pname(param_name) + if ( + f"{self._get_base_pname(param_name)}{int(param_values[0]):d}" + ) not in self.parameters: + raise ValueError( + f"Line {line_num:d} and {int(param_values[0]):d} : " + "Different reaction rules in parameter constraints." + ) + else: + if f"x[C.{param_name}" + ( + f"{line_num:d}]" if func_name != "user_defined" else "]" + ) != ( + f"x[C.{self._get_base_pname(param_name)}" + f"{int(param_values[0]):d}]" + ): + self.param_excluded.append( + f"{param_name}" + + (f"{line_num:d}" if func_name != "user_defined" else "") + ) + self.param_info.append( + f"x[C.{param_name}" + + (f"{line_num:d}]" if func_name != "user_defined" else "]") + + " = " + + f"x[C.{self._get_base_pname(param_name)}" + + f"{int(param_values[0]):d}]" + ) + self.param_constraints.append( + f"x[C.{param_name}" + + (f"{line_num:d}]" if func_name != "user_defined" else "]") + + " = " + + f"x[C.{self._get_base_pname(param_name)}" + + f"{int(param_values[0]):d}]" + ) + else: + raise ValueError( + f"line{line_num:d}: {line}\nInvalid expression in the input parameter." + ) + + def _process_ival_section(self, line_num: int, line: str) -> None: + + initial_values = line.split("|")[2].strip().split(",") + for ival in initial_values: + if ival.startswith("fixed "): + ival = ival.split("fixed ")[-1] + self.fixed_species.append(ival.split("=")[0].strip(" ")) + if ival.split("=")[0].strip(" ") in line.split("|")[0]: + if self._isfloat(ival.split("=")[1].strip(" ")): + self.init_info.append( + "y0[V." + + ival.split("=")[0].strip(" ") + + "] = " + + ival.split("=")[1].strip(" ") + ) + else: + raise ValueError(f"line{line_num:d}: Initial value must be int or float.") + else: + raise NameError( + f"line{line_num:d}: " f"Name'{ival.split('=')[0].strip(' ')}' is not defined." + ) + + @staticmethod + def _get_base_pname(param_name: str) -> str: + while param_name[-1].isdecimal(): + param_name = param_name[:-1] + return param_name + + def _preprocessing(self, func_name: str, line_num: int, line: str, *args: str) -> List[str]: + """ + Extract the information about parameter and/or initial values + if '|' in the line and find a keyword to identify reaction rules. + + Parameters + ---------- + func_name : str + Name of the rule function. + + line_num : int + Line number. + + line : str + Each line of the input text. + + Returns + ------- + description : list of strings + + """ + self._set_params(line_num, func_name, *args) + if "|" in line: + if line.split("|")[1].strip(): + self._process_pval_section(func_name, line_num, line, *args) + if line.count("|") > 1 and line.split("|")[2].strip(): + self._process_ival_section(line_num, line) + line = line.split("|")[0] + if func_name != "user_defined": + hit_words: List[str] = [] + for word in self.rule_words[func_name]: + # Choose longer word + if word in line: + hit_words.append(word) + description = line.strip().split(max(hit_words, key=len)) + if description[1] and not description[1].startswith(" "): + self._raise_detection_error(line_num, line) + else: + description = line.strip().split(":") + return description + + @staticmethod + def _word2scores(word: str, sentence: str) -> List[float]: + """ + Calculate similarity scores between word and sentence. + + Parameters + ---------- + word : str + User-defined word. + sentence : str + Textual unit consisting of two or more words. + + returns + ------- + scores : list + List containing similarity scores. + + """ + scores = [ + SequenceMatcher(None, word, sentence[i : i + len(word)]).ratio() + for i in range(len(sentence) - len(word) + 1) + ] + return scores + + def _get_partial_similarity(self, line: str) -> UnregisteredRule: + """ + Suggest similar rule word when user-defined word is not registered + in rule_words. + + Parameters + ---------- + line : str + Each line of the input text. + + Returns + ------- + unregistered_rule : UnregisteredRule + Rule word with the highest similarity score. + + """ + match_words = [] + match_scores = [] + str_subset = [] + for rules in self.rule_words.values(): + for word in rules: + ratio = self._word2scores(word, line) + if ratio: + match_words.append(word) + match_scores.append(max(ratio)) + str_subset.append(line[np.argmax(ratio) : np.argmax(ratio) + len(word)]) + expected_word = ( + None + if all([score < self.similarity_threshold for score in match_scores]) + else match_words[np.argmax(match_scores)] + ) + original_word = ( + None if expected_word is None else str_subset[match_words.index(expected_word)] + ) + unregistered_rule = UnregisteredRule(expected_word, original_word) + + return unregistered_rule + + @staticmethod + def _remove_prepositions(sentence: str) -> str: + """ + Remove preposition from text not to use it for identifying reaction rules. + """ + for preposition in PREPOSITIONS: + if sentence.endswith(f" {preposition}"): + return sentence[: -len(preposition) - 1] + return sentence + + def _get_arrow_error_message(self, line_num: int) -> str: + message = ( + f"line{line_num}: Use one of ({', '.join(self.fwd_arrows)}) for unidirectional " + f"reaction or ({', '.join(self.double_arrows)}) for bi-directional reaction" + ) + return message + + def _bind_and_dissociate(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'A + B --> AB' # bind, unidirectional + >>> 'AB --> A + B' # dissociate, unidirectional + >>> 'A + B <--> AB' # bind and dissociate, bidirectional + """ + for arrow in self._available_arrows(): + if arrow in line: + params_used = ["kf"] if arrow in self.fwd_arrows else ["kf", "kr"] + break + else: + raise ArrowError(self._get_arrow_error_message(line_num) + ".") + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, *params_used + ) + is_binding: bool + is_unidirectional: bool + for arrow in self._available_arrows(): + if arrow in description[1]: + is_binding = True + is_unidirectional = True if arrow in self.fwd_arrows else False + component1 = description[0].strip(" ") + component2 = description[1].split(arrow)[0].strip(" ") + complex = description[1].split(arrow)[1].strip(" ") + break + elif arrow in description[0]: + is_binding = False + is_unidirectional = True if arrow in self.fwd_arrows else False + component1 = description[0].split(arrow)[1].strip(" ") + component2 = description[1].strip(arrow) + complex = description[0].split(arrow)[0].strip(" ") + break + else: + raise ArrowError(self._get_arrow_error_message(line_num) + ".") + if component1 == complex or component2 == complex: + raise ValueError(f"line{line_num:d}: {complex} <- Use a different name.") + else: + self._set_species(component1, component2, complex) + self.complex_formations.append( + ComplexFormation(line_num, set([component1, component2]), complex, is_binding) + ) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{component1}] * y[V.{component2}]" + + (f" - x[C.kr{line_num:d}] * y[V.{complex}]" if not is_unidirectional else "") + if is_binding + else f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{complex}]" + + ( + f" - x[C.kr{line_num:d}] * y[V.{component1}] * y[V.{component2}]" + if not is_unidirectional + else "" + ) + ) + counter_component1, counter_component2, counter_complex = (0, 0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{component1}]" in eq: + counter_component1 += 1 + self.differential_equations[i] = ( + eq + (" - " if is_binding else " + ") + f"v[{line_num:d}]" + ) + elif f"dydt[V.{component2}]" in eq: + counter_component2 += 1 + self.differential_equations[i] = ( + eq + (" - " if is_binding else " + ") + f"v[{line_num:d}]" + ) + elif f"dydt[V.{complex}]" in eq: + counter_complex += 1 + self.differential_equations[i] = ( + eq + (" + " if is_binding else " - ") + f"v[{line_num:d}]" + ) + if counter_component1 == 0: + self.differential_equations.append( + f"dydt[V.{component1}] = - v[{line_num:d}]" + if is_binding + else f"dydt[V.{component1}] = + v[{line_num:d}]" + ) + if counter_component2 == 0: + self.differential_equations.append( + f"dydt[V.{component2}] = - v[{line_num:d}]" + if is_binding + else f"dydt[V.{component2}] = + v[{line_num:d}]" + ) + if counter_complex == 0: + self.differential_equations.append( + f"dydt[V.{complex}] = + v[{line_num:d}]" + if is_binding + else f"dydt[V.{complex}] = - v[{line_num:d}]" + ) + + def dimerize(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'A dimerizes <--> AA' + >>> 'A homodimerizes <--> AA' + >>> 'A forms a dimer <--> AA' + >>> 'A forms dimers <--> AA' + + Notes + ----- + * Parameters + .. math:: kf, kr + + * Rate equation + .. math:: v = kf * [A] * [A] - kr * [AA] + + * Differential equation + .. math:: + + d[A]]/dt = - 2 * v + + d[AA]/dt = + v + + """ + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, "kf", "kr" + ) + is_unidirectional: bool + monomer = description[0].strip(" ") + for arrow in self._available_arrows(): + if arrow in description[1]: + is_unidirectional = True if arrow in self.fwd_arrows else False + dimer = description[1].split(arrow)[1].strip(" ") + break + else: + raise ArrowError( + self._get_arrow_error_message(line_num) + " to specify the name of the dimer." + ) + if monomer == dimer: + raise ValueError(f"{dimer} <- Use a different name.") + self._set_species(monomer, dimer) + self.complex_formations.append(ComplexFormation(line_num, set(monomer), dimer, True)) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{monomer}] * y[V.{monomer}]" + + (f" - x[C.kr{line_num:d}] * y[V.{dimer}]" if not is_unidirectional else "") + ) + counter_monomer, counter_dimer = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{monomer}]" in eq: + counter_monomer += 1 + self.differential_equations[i] = eq + f" - 2 * v[{line_num:d}]" + elif f"dydt[V.{dimer}]" in eq: + counter_dimer += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_monomer == 0: + self.differential_equations.append(f"dydt[V.{monomer}] = - v[{line_num:d}]") + if counter_dimer == 0: + self.differential_equations.append(f"dydt[V.{dimer}] = + v[{line_num:d}]") + + def bind(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'A binds B <--> AB' + >>> 'A forms complexes with B <--> AB' + + Notes + ----- + * Parameters + .. math:: kf, kr + + * Rate equation + .. math:: v = kf * [A] * [B] - kr * [AB] + + * Differential equation + .. math:: + + d[A]/dt = - v + + d[B]/dt = - v + + d[AB]/dt = + v + + """ + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, "kf", "kr" + ) + is_unidirectional: bool + component1 = description[0].strip(" ") + for arrow in self._available_arrows(): + if arrow in description[1]: + is_unidirectional = True if arrow in self.fwd_arrows else False + component2 = description[1].split(arrow)[0].strip(" ") + complex = description[1].split(arrow)[1].strip(" ") + break + else: + raise ArrowError( + self._get_arrow_error_message(line_num) + + " to specify the name of the protein complex." + ) + if component1 == complex or component2 == complex: + raise ValueError(f"line{line_num:d}: {complex} <- Use a different name.") + elif component1 == component2: + self.dimerize(line_num, line) + else: + self._set_species(component1, component2, complex) + self.complex_formations.append( + ComplexFormation(line_num, set([component1, component2]), complex, True) + ) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{component1}] * y[V.{component2}]" + + (f" - x[C.kr{line_num:d}] * y[V.{complex}]" if not is_unidirectional else "") + ) + counter_component1, counter_component2, counter_complex = (0, 0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{component1}]" in eq: + counter_component1 += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{component2}]" in eq: + counter_component2 += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{complex}]" in eq: + counter_complex += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_component1 == 0: + self.differential_equations.append(f"dydt[V.{component1}] = - v[{line_num:d}]") + if counter_component2 == 0: + self.differential_equations.append(f"dydt[V.{component2}] = - v[{line_num:d}]") + if counter_complex == 0: + self.differential_equations.append(f"dydt[V.{complex}] = + v[{line_num:d}]") + + def dissociate(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'AB dissociates to A and B' + >>> 'AB is dissociated into A and B' + + Notes + ----- + * Parameters + .. math:: kf, kr + + * Rate equation + .. math:: v = kf * [AB] - kr * [A] * [B] + + * Differential equation + .. math:: + + d[A]/dt = + v + + d[B]/dt = + v + + d[AB]/dt = - v + + """ + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, "kf", "kr" + ) + complex = description[0].strip(" ") + if " and " not in description[1]: + raise ValueError( + f"Use 'and' in line{line_num:d}:\ne.g., AB is dissociated into A and B" + ) + else: + component1 = description[1].split(" and ")[0].strip(" ") + component2 = description[1].split(" and ")[1].strip(" ") + self._set_species(complex, component1, component2) + self.complex_formations.append( + ComplexFormation(line_num, set([component1, component2]), complex, False) + ) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{complex}]" + f" - x[C.kr{line_num:d}] * y[V.{component1}] * y[V.{component2}]" + ) + counter_complex, counter_component1, counter_component2 = (0, 0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{complex}]" in eq: + counter_complex += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{component1}]" in eq: + counter_component1 += 1 + self.differential_equations[i] = ( + eq + f" + v[{line_num:d}]" + if component1 != component2 + else eq + f" + 2 * v[{line_num:d}]" + ) + elif f"dydt[V.{component2}]" in eq: + counter_component2 += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_complex == 0: + self.differential_equations.append(f"dydt[V.{complex}] = - v[{line_num:d}]") + if counter_component1 == 0: + self.differential_equations.append(f"dydt[V.{component1}] = + v[{line_num:d}]") + if counter_component2 == 0: + self.differential_equations.append(f"dydt[V.{component2}] = + v[{line_num:d}]") + + def is_phosphorylated(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'uA is phosphorylated <--> pA' + + Notes + ----- + * Parameters + .. math:: kf, kr + + * Rate equation + .. math:: v = kf * [uA] - kr * [pA] + + * Differential equation + .. math:: + + d[uA]/dt = - v + + d[pA]/dt = + v + + """ + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, "kf", "kr" + ) + is_unidirectional: bool + unphosphorylated_form = description[0].strip(" ") + for arrow in self._available_arrows(): + if arrow in description[1]: + is_unidirectional = True if arrow in self.fwd_arrows else False + phosphorylated_form = description[1].split(arrow)[1].strip(" ") + break + else: + raise ArrowError( + self._get_arrow_error_message(line_num) + + " to specify the name of the phosphorylated protein." + ) + self._set_species(unphosphorylated_form, phosphorylated_form) + + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.kf{line_num:d}] * y[V.{unphosphorylated_form}]" + + ( + f" - x[C.kr{line_num:d}] * y[V.{phosphorylated_form}]" + if not is_unidirectional + else "" + ) + ) + counter_unphosphorylated_form, counter_phosphorylated_form = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{unphosphorylated_form}]" in eq: + counter_unphosphorylated_form += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif "dydt[V.{phosphorylated_form}]" in eq: + counter_phosphorylated_form += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_unphosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{unphosphorylated_form}] = - v[{line_num:d}]" + ) + if counter_phosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{phosphorylated_form}] = + v[{line_num:d}]" + ) + + def is_dephosphorylated(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'pA is dephosphorylated --> uA' + + Notes + ----- + * Parameters + .. math:: V, K + + * Rate equation + .. math:: v = V * [pA] / (K + [pA]) + + * Differential equation + .. math:: + + d[uA]/dt = + v + + d[pA]/dt = - v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "V", "K") + phosphorylated_form = description[0].strip(" ") + for arrow in self.fwd_arrows: + if arrow in description[1]: + unphosphorylated_form = description[1].split(arrow)[1].strip(" ") + break + else: + raise ArrowError( + f"line{line_num:d}: Use one of ({', '.join(self.fwd_arrows)}) " + "to specify the name of the dephosphorylated protein." + ) + self._set_species(phosphorylated_form, unphosphorylated_form) + + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.V{line_num:d}] * y[V.{phosphorylated_form}] / " + f"(x[C.K{line_num:d}] + y[V.{phosphorylated_form}])" + ) + counter_unphosphorylated_form, counter_phosphorylated_form = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{unphosphorylated_form}]" in eq: + counter_unphosphorylated_form += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + elif f"dydt[V.{phosphorylated_form}]" in eq: + counter_phosphorylated_form += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + if counter_unphosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{unphosphorylated_form}] = + v[{line_num:d}]" + ) + if counter_phosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{phosphorylated_form}] = - v[{line_num:d}]" + ) + + def phosphorylate(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'B phosphorylates uA --> pA' + + Notes + ----- + * Parameters + .. math:: V, K + + * Rate equation + .. math:: v = V * [B] * [uA] / (K + [uA]) + + * Differential equation + .. math:: + + d[uA]/dt = - v + + d[pA]/dt = + v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "V", "K") + kinase = description[0].strip(" ") + for arrow in self.fwd_arrows: + if arrow in description[1]: + unphosphorylated_form = description[1].split(arrow)[0].strip(" ") + phosphorylated_form = description[1].split(arrow)[1].strip(" ") + break + else: + raise ValueError( + f"line{line_num:d}: " + f"Use one of {', '.join(self.fwd_arrows)} to specify " + "the name of the phosphorylated (or activated) protein." + ) + if unphosphorylated_form == phosphorylated_form: + raise ValueError(f"line{line_num:d}: {phosphorylated_form} <- Use a different name.") + self._set_species(kinase, unphosphorylated_form, phosphorylated_form) + + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.V{line_num:d}] * y[V.{kinase}] * y[V.{unphosphorylated_form}] / " + f"(x[C.K{line_num:d}] + y[V.{unphosphorylated_form}])" + ) + counter_unphosphorylated_form, counter_phosphorylated_form = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{unphosphorylated_form}]" in eq: + counter_unphosphorylated_form += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{phosphorylated_form}]" in eq: + counter_phosphorylated_form += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_unphosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{unphosphorylated_form}] = - v[{line_num:d}]" + ) + if counter_phosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{phosphorylated_form}] = + v[{line_num:d}]" + ) + + def dephosphorylate(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'B dephosphorylates pA --> uA' + + Notes + ----- + * Parameters + .. math:: V, K + + * Rate equation + .. math:: v = V * [B] * [pA] / (K + [pA]) + + * Differential equation + .. math:: + + d[uA]/dt = + v + + d[pA]/dt = - v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "V", "K") + phosphatase = description[0].strip(" ") + for arrow in self.fwd_arrows: + if arrow in description[1]: + phosphorylated_form = description[1].split(arrow)[0].strip(" ") + unphosphorylated_form = description[1].split(arrow)[1].strip(" ") + break + else: + raise ArrowError( + f"line{line_num:d}: " + f"Use one of {', '.join(self.fwd_arrows)} to specify " + "the name of the dephosphorylated (or deactivated) protein." + ) + if phosphorylated_form == unphosphorylated_form: + raise ValueError(f"line{line_num:d}: {unphosphorylated_form} <- Use a different name.") + self._set_species(phosphatase, phosphorylated_form, unphosphorylated_form) + + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.V{line_num:d}] * y[V.{phosphatase}] * y[V.{phosphorylated_form}] / " + f"(x[C.K{line_num:d}] + y[V.{phosphorylated_form}])" + ) + counter_phosphorylated_form, counter_unphosphorylated_form = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{phosphorylated_form}]" in eq: + counter_phosphorylated_form += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{unphosphorylated_form}]" in eq: + counter_unphosphorylated_form += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_phosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{phosphorylated_form}] = - v[{line_num:d}]" + ) + if counter_unphosphorylated_form == 0: + self.differential_equations.append( + f"dydt[V.{unphosphorylated_form}] = + v[{line_num:d}]" + ) + + def transcribe(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'B transcribes a' + >>> 'B1 & B2 transcribe a' # (AND-gate) + >>> 'B transcribes a, repressed by C' # (Negative regulation) + + Notes + ----- + * Parameters + .. math:: V, K, n, (KF, nF) + + * Rate equation + .. math:: + + v = V * [B] ^ {n} / (K ^ {n} + [B] ^ {n}) + + v = V * ([B1] * [B2]) ^ {n} / (K ^ {n} + ([B1] * [B2]) ^ {n}) + + v = V * [B] ^ {n} / (K ^ {n} + [B] ^ {n} + ([C] / KF) ^ {nF}) + + * Differential equation + .. math:: d[a]/dt = + v + + """ + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, "V", "K", "n", "KF", "nF" + ) + repressor: Optional[str] = None + ratio = self._word2scores(", repressed by", description[1]) + if not ratio or max(ratio) < 1.0: + self.parameters.remove(f"KF{line_num:d}") + self.parameters.remove(f"nF{line_num:d}") + mRNA = description[1].strip() + if " " in mRNA: + # Fix typo in line{line_num:d} + raise ValueError( + f"line{line_num:d}: " + "Add ', repressed by XXX' to describe negative regulation from XXX." + ) + else: + # Add negative regulation from repressor + mRNA = description[1].split(", repressed by")[0].strip() + repressor = description[1].split(", repressed by")[1].strip() + if " & " not in description[0]: + TF = description[0].strip(" ") + self._set_species(mRNA, TF) + if repressor is not None: + self._set_species(repressor) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.V{line_num:d}] * y[V.{TF}] ** x[C.n{line_num:d}] / " + f"(x[C.K{line_num:d}] ** x[C.n{line_num:d}] + " + f"y[V.{TF}] ** x[C.n{line_num:d}]" + + ( + ")" + if repressor is None + else f" + (y[V.{repressor}] / x[C.KF{line_num:d}]) ** x[C.nF{line_num:d}])" + ) + ) + else: + # AND-gate + TFs = [TF.strip(" ") for TF in description[0].split(" & ")] + self._set_species(mRNA, *TFs) + if repressor is not None: + self._set_species(repressor) + self.reactions.append( + f"v[{line_num:d}] = " + f"x[C.V{line_num:d}] * ({'y[V.' + '] * y[V.'.join(TFs) + ']'}) ** x[C.n{line_num:d}] / " + f"(x[C.K{line_num:d}] ** x[C.n{line_num:d}] + " + f"({'y[V.' + '] * y[V.'.join(TFs) + ']'}) ** x[C.n{line_num:d}]" + + ( + ")" + if repressor is None + else f" + (y[V.{repressor}] / x[C.KF{line_num:d}]) ** x[C.nF{line_num:d}])" + ) + ) + counter_mRNA = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{mRNA}]" in eq: + counter_mRNA += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_mRNA == 0: + self.differential_equations.append(f"dydt[V.{mRNA}] = + v[{line_num:d}]") + + def synthesize(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'B synthesizes A' + + Notes + ----- + * Parameters + .. math:: kf + + * Rate equation + .. math:: v = kf * [B] + + * Differential equation + .. math:: d[A]/dt = + v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "kf") + catalyst = description[0].strip(" ") + product = description[1].strip(" ") + self._set_species(catalyst, product) + self.reactions.append(f"v[{line_num:d}] = x[C.kf{line_num:d}] * y[V.{catalyst}]") + counter_product = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{product}]" in eq: + counter_product += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_product == 0: + self.differential_equations.append(f"dydt[V.{product}] = + v[{line_num:d}]") + + def is_synthesized(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'A is synthesized' + + Notes + ----- + * Parameters + .. math:: kf + + * Rate equation + .. math:: v = kf + + * Differential equation + .. math:: d[A]/dt = + v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "kf") + chemical_species = description[0].strip(" ") + self._set_species(chemical_species) + self.reactions.append(f"v[{line_num:d}] = x[C.kf{line_num:d}]") + counter_chemical_species = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{chemical_species}]" in eq: + counter_chemical_species += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_chemical_species == 0: + self.differential_equations.append(f"dydt[V.{chemical_species}] = + v[{line_num:d}]") + + def degrade(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'B degrades A' + + Notes + ----- + * Parameters + .. math:: kf + + * Rate equation + .. math:: v = kf * [B] * [A] + + * Differential equation + .. math:: d[A]/dt = - v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "kf") + protease = description[0].strip(" ") + protein = description[1].strip(" ") + self._set_species(protease, protein) + self.reactions.append( + f"v[{line_num:d}] = x[C.kf{line_num:d}] * y[V.{protease}] * y[V.{protein}]" + ) + counter_protein = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{protein}]" in eq: + counter_protein += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + if counter_protein == 0: + self.differential_equations.append(f"dydt[V.{protein}] = - v[{line_num:d}]") + + def is_degraded(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'A is degraded' + + Notes + ----- + * Parameters + .. math:: kf + + * Rate equation + .. math:: v = kf * [A] + + * Differential equation + .. math:: d[A]/dt = - v + + """ + description = self._preprocessing(sys._getframe().f_code.co_name, line_num, line, "kf") + chemical_species = description[0].strip(" ") + self._set_species(chemical_species) + self.reactions.append(f"v[{line_num:d}] = x[C.kf{line_num:d}] * y[V.{chemical_species}]") + counter_chemical_species = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{chemical_species}]" in eq: + counter_chemical_species += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + if counter_chemical_species == 0: + self.differential_equations.append(f"dydt[V.{chemical_species}] = - v[{line_num:d}]") + + def translocate(self, line_num: int, line: str) -> None: + r""" + Examples + -------- + >>> 'A_at_cyt translocates from cytoplasm to nucleus (V_cyt, V_nuc) <--> A_at_nuc' + >>> 'A_at_cyt is translocated from cytoplasm to nucleus (V_cyt, V_nuc) <--> A_at_nuc' + + Notes + ----- + * Parameters + .. math:: kf, kr, (V_{pre}, V_{post}) + + * Rate equation + .. math:: v = kf * [A\_at\_pre] - kr * (V_{post} / V_{pre}) * [A\_at\_post] + + * Differential equation + .. math:: + + d[A\_at\_pre]/dt = - v + + d[A\_at\_post]/dt = + v * (V_{pre} / V_{post}) + + """ + for arrow in self._available_arrows(): + if arrow in line: + is_unidirectional = True if arrow in self.fwd_arrows else False + params_used = ["kf"] if is_unidirectional else ["kf", "kr"] + break + else: + raise ArrowError(self._get_arrow_error_message(line_num)) + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, *params_used + ) + pre_translocation = description[0].strip(" ") + for arrow in self._available_arrows(): + if arrow in description[1]: + post_translocation = description[1].split(arrow)[1].strip(" ") + break + else: + assert False + if pre_translocation == post_translocation: + raise ValueError(f"line{line_num:d}: {post_translocation} <- Use a different name.") + # Information about compartment volumes + if "(" in description[1] and ")" in description[1]: + [pre_volume, post_volume] = description[1].split("(")[-1].split(")")[0].split(",") + if not self._isfloat(pre_volume.strip(" ")) or not self._isfloat( + post_volume.strip(" ") + ): + raise ValueError("pre_volume and post_volume must be float or int.") + else: + [pre_volume, post_volume] = ["1", "1"] + self._set_species(pre_translocation, post_translocation) + self.reactions.append( + f"v[{line_num:d}] = x[C.kf{line_num:d}] * y[V.{pre_translocation}]" + + ( + f" - x[C.kr{line_num:d}] * y[V.{post_translocation}]" + if not is_unidirectional + else "" + ) + ) + if float(pre_volume.strip(" ")) != float(post_volume.strip(" ")): + self.reactions[ + -1 + ] = f"v[{line_num:d}] = " f"x[C.kf{line_num:d}] * y[V.{pre_translocation}]" + ( + f" - x[C.kr{line_num:d}] * " + f"({post_volume.strip()} / {pre_volume.strip()}) * " + f"y[V.{post_translocation}]" + if not is_unidirectional + else "" + ) + counter_pre_translocation, counter_post_translocation = (0, 0) + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{pre_translocation}]" in eq: + counter_pre_translocation += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{post_translocation}]" in eq: + counter_post_translocation += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if float(pre_volume.strip(" ")) != float(post_volume.strip(" ")): + self.differential_equations[ + i + ] += f" * ({pre_volume.strip()} / {post_volume.strip()})" + if counter_pre_translocation == 0: + self.differential_equations.append(f"dydt[V.{pre_translocation}] = - v[{line_num:d}]") + if counter_post_translocation == 0: + self.differential_equations.append(f"dydt[V.{post_translocation}] = + v[{line_num:d}]") + if float(pre_volume.strip(" ")) != float(post_volume.strip(" ")): + self.differential_equations[ + -1 + ] += f" * ({pre_volume.strip()} / {post_volume.strip()})" + + def user_defined(self, line_num: int, line: str) -> None: + """ + Examples + -------- + >>> 'Reactant --> Product: define rate equation here' + + Notes + ----- + * Use p[xxx] and u[xxx] for describing parameters and species, respectively. + + * Use '0' or '∅' for degradation/creation to/from nothing. + + * Differential equation + .. math:: + + d[Reactant]/dt = - v + + d[Product]/dt = + v + """ + all_params = re.findall(r"p\[(.*?)\]", line) + description = self._preprocessing( + sys._getframe().f_code.co_name, line_num, line, *all_params + ) + balance = description[0].strip() + rate_equation = description[1].strip() + for arrow in self.fwd_arrows: + if arrow in balance: + reactant, product = balance.split(arrow) + self._set_species(reactant.strip(), product.strip()) + break + else: + raise ArrowError(f"line{line_num:d}: Use one of {', '.join(self.fwd_arrows)}.") + rate_equation = ( + rate_equation.replace("p[", "x[C.").replace("u[", "y[V.").replace("^", "**") + ) + self.reactions.append(f"v[{line_num:d}] = " + rate_equation.strip()) + counter_reactant = 0 + counter_product = 0 + for i, eq in enumerate(self.differential_equations): + if f"dydt[V.{reactant}]" in eq and reactant not in self.nothing: + counter_reactant += 1 + self.differential_equations[i] = eq + f" - v[{line_num:d}]" + elif f"dydt[V.{product}]" in eq and product not in self.nothing: + counter_product += 1 + self.differential_equations[i] = eq + f" + v[{line_num:d}]" + if counter_reactant == 0 and reactant not in self.nothing: + self.differential_equations.append(f"dydt[V.{reactant}] = - v[{line_num:d}]") + if counter_product == 0 and product not in self.nothing: + self.differential_equations.append(f"dydt[V.{product}] = + v[{line_num:d}]") + + def _extract_event(self, line_num: int, line: str): + # About biochemical event + if line.startswith("@rxn "): + line = self._remove_prefix(line, "@rxn ") + if line.count(":") != 1: + raise SyntaxError(f"line{line_num:d}: Missing colon") + else: + self.user_defined(line_num, line) + # About observables + elif line.startswith("@obs "): + line = self._remove_prefix(line, "@obs ") + if line.count(":") != 1: + raise SyntaxError( + f"line{line_num:d}: Missing colon\n" + "Should be `@obs : `." + ) + else: + self.obs_desc.append(line.split(":")) + # About simulation info. + elif line.startswith("@sim "): + line = self._remove_prefix(line, "@sim ") + if line.count(":") != 1: + raise SyntaxError(f"line{line_num:d}: Missing colon") + else: + if line.startswith("tspan"): + t_info = line.split(":")[-1].strip() + if "[" in t_info and "]" in t_info: + [t0, tf] = t_info.split("[")[-1].split("]")[0].split(",") + if t0.strip(" ").isdecimal() and tf.strip(" ").isdecimal(): + self.sim_tspan.append(t0) + self.sim_tspan.append(tf) + else: + raise TypeError("@sim tspan: [t0, tf] must be a list of integers.") + else: + raise ValueError( + "`tspan` must be a two element vector [t0, tf] " + "specifying the initial and final times." + ) + elif line.startswith("unperturbed"): + self.sim_unperturbed += line.split(":")[-1].strip() + elif line.startswith("condition "): + self.sim_conditions.append(self._remove_prefix(line, "condition ").split(":")) + else: + raise ValueError( + f"(line{line_num:d}) Available options are: " + "'@sim tspan:', '@sim unperturbed:', or '@sim condition XXX:'." + ) + # Additional species + elif line.startswith("@add "): + line = self._remove_prefix(line, "@add ") + if line.startswith("species "): + line = self._remove_prefix(line, "species ") + new_species = line.strip() + if new_species not in self.species: + self._set_species(new_species) + else: + raise NameError(f"{new_species} is already defined.") + elif line.startswith("param "): + line = self._remove_prefix(line, "param ") + new_param = line.strip() + if new_param not in self.parameters: + self._set_params(None, None, new_param) + self.param_excluded.append(new_param) + else: + raise NameError(f"{new_param} is already defined.") + else: + raise ValueError(f"(line{line_num:d}) Must be either @add param or @add species.") + else: + raise ValueError("Available symbols are: @rxn, @add, @obs, @sim.") + + def create_ode(self) -> None: + """ + Find a keyword in each line to identify the reaction rule and + construct an ODE model. + + """ + with open(self.input_txt, encoding="utf-8") as f: + lines = f.readlines() + for line_num, line in enumerate(lines, start=1): + # Remove double spaces + while True: + if " " not in line: + break + else: + line = line.replace(" ", " ") + # Comment out + line = line.split("#")[0].rstrip(" ") + if not line.strip(): + # Skip blank lines + continue + elif lines.count(line) > 1: + # Find duplicate lines + raise DuplicateError( + f"Reaction '{line}' is duplicated in lines " + + ", ".join([str(i + 1) for i, rxn in enumerate(lines) if rxn == line]) + ) + elif line.startswith("@"): + self._extract_event(line_num, line) + # Detect reaction rule + else: + for reaction_rule, words in self.rule_words.items(): + if any([self._remove_prepositions(word) in line for word in words]): + exec("self." + reaction_rule + "(line_num, line)") + break + else: + self._raise_detection_error(line_num, line) diff --git a/biomass/construction/template/README.md b/biomass/construction/template/README.md new file mode 100644 index 00000000..732c1225 --- /dev/null +++ b/biomass/construction/template/README.md @@ -0,0 +1,15 @@ +# Template for BioMASS model construction + +A brief description of each file/folder is below: + +| Name | Content | +| ---------------------------------------------- | -------------------------------------------------------------------------------------------------------- | +| [`name2idx/`](./name2idx/) | Names of model parameters and species | +| [`reaction_network.py`](./reaction_network.py) | Flux vector and reaction indices grouped according to biological processes | +| [`ode.py`](./ode.py) | Differential equation, parameters and initial condition | +| [`observalbe.py`](./observable.py) | Observables, simulations and experimental data | +| [`viz.py`](./viz.py) | Plotting parameters for customizing figure properties | +| [`search_param.py`](./search_param.py) | Lower and upper bounds of model parameters to be estimated | +| [`problem.py`](./problem.py) | An objective function to be minimized, i.e., the distance between model simulation and experimental data | + +Example models can be found at https://biomass-core.readthedocs.io/en/latest/models.html. \ No newline at end of file diff --git a/biomass/construction/template/__init__.py b/biomass/construction/template/__init__.py new file mode 100644 index 00000000..a0bc1633 --- /dev/null +++ b/biomass/construction/template/__init__.py @@ -0,0 +1,5 @@ +from .name2idx import C, V +from .ode import initial_values, param_values +from .problem import OptimizationProblem +from .reaction_network import ReactionNetwork +from .viz import Visualization diff --git a/biomass/construction/template/name2idx/__init__.py b/biomass/construction/template/name2idx/__init__.py new file mode 100644 index 00000000..3869ddc0 --- /dev/null +++ b/biomass/construction/template/name2idx/__init__.py @@ -0,0 +1,2 @@ +from .parameters import C +from .species import V diff --git a/biomass/construction/template/name2idx/parameters.py b/biomass/construction/template/name2idx/parameters.py new file mode 100644 index 00000000..cb3917c0 --- /dev/null +++ b/biomass/construction/template/name2idx/parameters.py @@ -0,0 +1,19 @@ +from dataclasses import make_dataclass +from typing import Dict, List + +NAMES: List[str] = [] + +NUM: int = len(NAMES) + +Parameters = make_dataclass( + cls_name="Parameters", + fields=[(name, int) for name in NAMES], + namespace={"NAMES": NAMES, "NUM": NUM}, + frozen=True, +) + +name2idx: Dict[str, int] = {k: v for v, k in enumerate(NAMES)} + +C = Parameters(**name2idx) + +del name2idx diff --git a/biomass/construction/template/name2idx/species.py b/biomass/construction/template/name2idx/species.py new file mode 100644 index 00000000..81096706 --- /dev/null +++ b/biomass/construction/template/name2idx/species.py @@ -0,0 +1,19 @@ +from dataclasses import make_dataclass +from typing import Dict, List + +NAMES: List[str] = [] + +NUM: int = len(NAMES) + +Species = make_dataclass( + cls_name="Species", + fields=[(name, int) for name in NAMES], + namespace={"NAMES": NAMES, "NUM": NUM}, + frozen=True, +) + +name2idx: Dict[str, int] = {k: v for v, k in enumerate(NAMES)} + +V = Species(**name2idx) + +del name2idx diff --git a/biomass/construction/template/observable.py b/biomass/construction/template/observable.py new file mode 100644 index 00000000..f990c233 --- /dev/null +++ b/biomass/construction/template/observable.py @@ -0,0 +1,78 @@ +from typing import List + +import numpy as np + +from biomass.dynamics.solver import * + +from .name2idx import C, V +from .ode import DifferentialEquation + + +class Observable(DifferentialEquation): + """ + Correlating model simulations and experimental measurements. + + Attributes + ---------- + obs_names : list of strings + Names of model observables. + + t : range + Simulation time span. + + conditions : list of strings + Experimental conditions. + + simulations : numpy.ndarray + The numpy array to store simulation results. + + normalization : nested dict + * 'timepoint' : Optional[int] + The time point at which simulated values are normalized. + If :obj:`None`, the maximum value will be used for normalization. + + * 'condition' : list of strings + The experimental conditions to use for normalization. + If empty, all conditions defined in ``sim.conditions`` will be used. + + experiments : list of dict + Time series data. + + error_bars : list of dict + Error bars to show in figures. + + """ + + def __init__(self): + super(Observable, self).__init__(perturbation={}) + self.obs_names: list = [] + self.t: range = range(101) + self.conditions: list = [] + self.simulations: np.ndarray = np.empty( + (len(self.obs_names), len(self.conditions), len(self.t)) + ) + self.normalization: dict = {} + self.experiments: list = [None] * len(self.obs_names) + self.error_bars: list = [None] * len(self.obs_names) + + def simulate(self, x: list, y0: list, _perturbation: dict = {}): + if _perturbation: + self.perturbation = _perturbation + # unperturbed steady state + + for i, condition in enumerate(self.conditions): + + sol = solve_ode(self.diffeq, y0, self.t, tuple(x)) + + if sol is None: + return False + else: + pass + + def set_data(self) -> None: + pass + + def get_timepoint(self, obs_name: str) -> List[int]: + if obs_name in self.obs_names: + return [] + assert False diff --git a/biomass/construction/template/ode.py b/biomass/construction/template/ode.py new file mode 100644 index 00000000..3475633e --- /dev/null +++ b/biomass/construction/template/ode.py @@ -0,0 +1,34 @@ +from .name2idx import C, V +from .reaction_network import ReactionNetwork + + +class DifferentialEquation(ReactionNetwork): + def __init__(self, perturbation): + super(DifferentialEquation, self).__init__() + self.perturbation = perturbation + + def diffeq(self, t, y, *x): + """Kinetic equations""" + v = self.flux(t, y, x) + + if self.perturbation: + for i, dv in self.perturbation.items(): + v[i] = v[i] * dv + + dydt = [0] * V.NUM + + return dydt + + +def param_values(): + """Parameter values""" + x = [1] * C.NUM + + return x + + +def initial_values(): + """Values of the initial condition""" + y0 = [0] * V.NUM + + return y0 diff --git a/biomass/construction/template/problem.py b/biomass/construction/template/problem.py new file mode 100644 index 00000000..90ef8635 --- /dev/null +++ b/biomass/construction/template/problem.py @@ -0,0 +1,103 @@ +import numpy as np +from scipy.spatial.distance import cosine + +from .observable import Observable +from .search_param import SearchParam + + +class OptimizationProblem(Observable, SearchParam): + def __init__(self): + super(OptimizationProblem, self).__init__() + + @property + def bounds(self): + """ + Lower and upper bounds on independent variables. + """ + search_region: np.ndarray = self.get_region() + lb = 10 ** search_region[0] + ub = 10 ** search_region[1] + return tuple(zip(lb, ub)) + + @staticmethod + def _compute_objval_rss(sim_data, exp_data): + """Return Residual Sum of Squares""" + return np.dot((sim_data - exp_data), (sim_data - exp_data)) + + @staticmethod + def _compute_objval_cos(sim_data, exp_data): + """Return Cosine distance""" + return cosine(sim_data, exp_data) + + @staticmethod + def _diff_sim_and_exp(sim_matrix, exp_dict, exp_timepoint, conditions, sim_norm_max): + sim_val = [] + exp_val = [] + + for idx, condition in enumerate(conditions): + if condition in exp_dict.keys(): + sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint))]) + exp_val.extend(exp_dict[condition]) + + return np.array(sim_val) / sim_norm_max, np.array(exp_val) + + def objective(self, indiv, *args): + """Define an objective function to be minimized.""" + if len(args) == 0: + (x, y0) = self.update(indiv) + elif len(args) == 1: + raise ValueError("not enough values to unpack (expected 2, got 1)") + elif len(args) == 2: + (x, y0) = args + else: + raise ValueError("too many values to unpack (expected 2)") + + self.set_data() + + if self.simulate(x, y0) is None: + error = np.zeros(len(self.obs_names)) + for i, obs_name in enumerate(self.obs_names): + if self.experiments[i] is not None: + error[i] = self._compute_objval_rss( + *self._diff_sim_and_exp( + self.simulations[i], + self.experiments[i], + self.get_timepoint(obs_name), + self.conditions, + sim_norm_max=1 + if not self.normalization + else ( + np.max( + self.simulations[ + self.obs_names.index(obs_name), + [ + self.conditions.index(c) + for c in ( + self.normalization[obs_name]["condition"] + if self.normalization[obs_name]["condition"] + else self.conditions + ) + ], + self.normalization[obs_name]["timepoint"], + ] + ) + if self.normalization[obs_name]["timepoint"] is not None + else np.max( + self.simulations[ + self.obs_names.index(obs_name), + [ + self.conditions.index(c) + for c in ( + self.normalization[obs_name]["condition"] + if self.normalization[obs_name]["condition"] + else self.conditions + ) + ], + ] + ) + ), + ) + ) + return np.sum(error) # < 1e12 + else: + return 1e12 diff --git a/biomass/construction/template/reaction_network.py b/biomass/construction/template/reaction_network.py new file mode 100644 index 00000000..f5523ed3 --- /dev/null +++ b/biomass/construction/template/reaction_network.py @@ -0,0 +1,39 @@ +from typing import Dict, List + +from biomass.analysis.reaction import is_duplicate + +from .name2idx import C, V + + +class ReactionNetwork(object): + """ + Reaction indices grouped according to biological processes. + This is used for sensitivity analysis (``target``='reaction'). + """ + + def __init__(self) -> None: + super(ReactionNetwork, self).__init__() + self.reactions: Dict[str, List[int]] = {} + + def group(self): + """ + Group reactions according to biological processes. + """ + biological_processes = [] + for process, indices in self.reactions.items(): + if not isinstance(indices, list): + raise TypeError("Use list for reaction indices in {}".format(process)) + biological_processes.append(indices) + + if not is_duplicate(self.reactions, biological_processes): + return biological_processes + + @staticmethod + def flux(t, y, x): + """ + Flux vector. + """ + + v = {} + + return v diff --git a/biomass/construction/template/search_param.py b/biomass/construction/template/search_param.py new file mode 100644 index 00000000..602d6a67 --- /dev/null +++ b/biomass/construction/template/search_param.py @@ -0,0 +1,78 @@ +import numpy as np + +from biomass.estimation import convert_scale, initialize_search_param + +from .name2idx import C, V +from .ode import initial_values, param_values + + +class SearchParam(object): + """Specify model parameters and/or initial values to optimize.""" + + def __init__(self): + # parameters + self.idx_params = [] + + # initial values + self.idx_initials = [] + + def get_region(self): + x = param_values() + y0 = initial_values() + + search_param = initialize_search_param( + parameters=C.NAMES, + species=V.NAMES, + param_values=x, + initial_values=y0, + estimated_params=self.idx_params, + estimated_initials=self.idx_initials, + ) + + search_rgn = np.zeros((2, len(x) + len(y0))) + # Default: 0.1 ~ 10 + for i, j in enumerate(self.idx_params): + search_rgn[0, j] = search_param[i] * 0.1 # lower bound + search_rgn[1, j] = search_param[i] * 10.0 # upper bound + # Default: 0.5 ~ 2 + for i, j in enumerate(self.idx_initials): + search_rgn[0, j + len(x)] = search_param[i + len(self.idx_params)] * 0.5 # lower bound + search_rgn[1, j + len(x)] = search_param[i + len(self.idx_params)] * 2.0 # upper bound + + # search_rgn[:,C.parameter] = [lower_bound, upper_bound] + # search_rgn[:,V.specie+len(x)] = [lower_bound, upper_bound] + + search_rgn = convert_scale( + region=search_rgn, + parameters=C.NAMES, + species=V.NAMES, + estimated_params=self.idx_params, + estimated_initials=self.idx_initials, + ) + + return search_rgn + + def update(self, indiv): + x = param_values() + y0 = initial_values() + + for i, j in enumerate(self.idx_params): + x[j] = indiv[i] + for i, j in enumerate(self.idx_initials): + y0[j] = indiv[i + len(self.idx_params)] + + # parameter constraints + + return x, y0 + + def gene2val(self, indiv_gene): + search_rgn = self.get_region() + indiv = 10 ** (indiv_gene * (search_rgn[1, :] - search_rgn[0, :]) + search_rgn[0, :]) + + return indiv + + def val2gene(self, indiv): + search_rgn = self.get_region() + indiv_gene = (np.log10(indiv) - search_rgn[0, :]) / (search_rgn[1, :] - search_rgn[0, :]) + + return indiv_gene diff --git a/biomass/construction/template/viz.py b/biomass/construction/template/viz.py new file mode 100644 index 00000000..eb6026bf --- /dev/null +++ b/biomass/construction/template/viz.py @@ -0,0 +1,79 @@ +from typing import List + +from matplotlib import pyplot as plt + +from biomass.plotting import * + +from .observable import Observable + + +class Visualization(Observable): + """ + Plotting parameters for customizing figure properties. + + Attributes + ---------- + cm : matplotlib.colors.ListedColormap (default: ``plt.cm.get_cmap('tab10')``) + Choosing colormaps for ``cmap``. + single_observable_options : list of SingleObservable + Visualization options for time-course simulation (single-observable). + multiple_observables_options : MultipleObservables + Visualization options for time-course simulation (multi-observables). + sensitivity_options : SensitivityOptions + Visualization options for sensitivity analysis results. + """ + + def __init__(self) -> None: + super().__init__() + + self.cm = plt.cm.get_cmap("tab10") + self.single_observable_options = [ + SingleObservable(self.cm, obs_name) for obs_name in self.obs_names + ] + self.multiple_observables_options = MultipleObservables(self.cm) + self.sensitivity_options = SensitivityOptions(self.cm) + + def get_single_observable_options(self) -> List[SingleObservable]: + + return self.single_observable_options + + def get_multiple_observables_options(self) -> MultipleObservables: + + return self.multiple_observables_options + + def get_sensitivity_options(self) -> SensitivityOptions: + + return self.sensitivity_options + + @staticmethod + def set_timecourse_rcParams() -> None: + """figure/simulation""" + plt.rcParams["font.size"] = 12 + plt.rcParams["axes.linewidth"] = 1.5 + plt.rcParams["xtick.major.width"] = 1.5 + plt.rcParams["ytick.major.width"] = 1.5 + plt.rcParams["lines.linewidth"] = 1.8 + plt.rcParams["lines.markersize"] = 12 + plt.rcParams["savefig.bbox"] = "tight" + # plt.rcParams["savefig.format"] = "pdf" + # plt.rcParams['font.family'] = 'Arial' + # plt.rcParams['mathtext.fontset'] = 'custom' + # plt.rcParams['mathtext.it'] = 'Arial:italic' + + @staticmethod + def set_sensitivity_rcParams() -> None: + """figure/sensitivity""" + plt.rcParams["font.size"] = 12 + plt.rcParams["axes.linewidth"] = 1.2 + plt.rcParams["xtick.major.width"] = 1.2 + plt.rcParams["ytick.major.width"] = 1.2 + plt.rcParams["savefig.bbox"] = "tight" + # plt.rcParams["savefig.format"] = "pdf" + # plt.rcParams['font.family'] = 'Arial' + + @staticmethod + def convert_species_name(name: str) -> str: + """figure/sensitivity/initial_condition + - Sensitivity for species with nonzero initial conditions. + """ + return name diff --git a/biomass/construction/text2model.py b/biomass/construction/text2model.py new file mode 100644 index 00000000..cf287ff3 --- /dev/null +++ b/biomass/construction/text2model.py @@ -0,0 +1,904 @@ +import os +import re +import shutil +from dataclasses import dataclass, field +from typing import Dict, Final, List, Literal, Optional + +from . import julia_template as jl +from .reaction_rules import ReactionRules + + +@dataclass +class Text2Model(ReactionRules): + """ + Build a BioMASS-formatted model based on template. + + **reaction** | **parameters** | **initial conditions** + + Attributes + ---------- + input_txt : str + Model description file (*.txt), e.g., 'Kholodenko1999.txt' + similarity_threshold : float (default: 0.7) + Similarity threshold used in text-to-model conversion. Must lie within (0, 1). + lang : Literal["python", "julia"] (default: 'python') + Either 'python' or 'julia'. + + - 'python': biomass (https://github.com/biomass-dev/biomass) + - 'julia': BioMASS.jl (https://github.com/biomass-dev/BioMASS.jl) + """ + + input_txt: str + similarity_threshold: float = 0.7 + lang: Literal["python", "julia"] = "python" + indentation: Final[str] = field(default=4 * " ", init=False) + + def __post_init__(self) -> None: + if not os.path.isfile(self.input_txt): + raise FileNotFoundError(f"{self.input_txt} does not exist.") + if self.lang not in ["python", "julia"]: + raise ValueError("lang must be either 'python' or 'julia'.") + self.name: str = os.path.splitext(self.input_txt)[0] + + def _update_parameters(self) -> None: + """ + Update name2idx/parameters.py + """ + if self.lang == "python": + with open( + os.path.join( + os.path.dirname(__file__), + "template", + "name2idx", + "parameters.py", + ), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + for line_num, line in enumerate(lines): + if line.startswith("NAMES: List[str] = []"): + lines[line_num] = "NAMES: List[str] = [\n" + lines[line_num] += ( + f'{self.indentation}"' + + f'",\n{self.indentation}"'.join(self.parameters) + + '",\n' + ) + lines[line_num] += "]\n" + with open( + os.path.join( + f"{self.name}", + "name2idx", + "parameters.py", + ), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + else: + lines = jl.PARAMETERS.splitlines() + for line_num, line in enumerate(lines): + if line.startswith("const NAMES = []"): + lines[line_num] = "const NAMES = [\n" + lines[line_num] += ( + f'{self.indentation}"' + + f'",\n{self.indentation}"'.join(self.parameters) + + '",\n' + ) + lines[line_num] += "]\n" + with open( + os.path.join( + f"{self.name}_jl", + "name2idx", + "parameters.jl", + ), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + + def _update_species(self) -> None: + """ + Update name2idx/species.py + """ + if self.lang == "python": + with open( + os.path.join( + os.path.dirname(__file__), + "template", + "name2idx", + "species.py", + ), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + for line_num, line in enumerate(lines): + if line.startswith("NAMES: List[str] = []"): + lines[line_num] = "NAMES: List[str] = [\n" + lines[line_num] += ( + f'{self.indentation}"' + + '",\n{}"'.format(self.indentation).join(self.species) + + '",\n' + ) + lines[line_num] += "]\n" + with open( + os.path.join( + f"{self.name}", + "name2idx", + "species.py", + ), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + else: + lines = jl.SPECIES.splitlines() + for line_num, line in enumerate(lines): + if line.startswith("const NAMES = []"): + lines[line_num] = "const NAMES = [\n" + lines[line_num] += ( + '{}"'.format(self.indentation) + + '",\n{}"'.format(self.indentation).join(self.species) + + '",\n' + ) + lines[line_num] += "]\n" + with open( + os.path.join( + f"{self.name}_jl", + "name2idx", + "species.jl", + ), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + + def _update_diffeq(self) -> None: + """Update differential equations. + + * Set flux vector (v). + * Set rhs of ODE. + + If you don't write info., all parameter values and initial values are + initialized to 1.0 and 0.0, respectively. + """ + if self.lang == "python": + with open( + os.path.join(os.path.dirname(__file__), "template", "ode.py"), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + # Check species which should be held fixed during simulation. + for species in self.fixed_species: + for i, equation in enumerate(self.differential_equations): + if f"dydt[V.{species}] = " in equation: + self.differential_equations[i] = equation.replace( + f"dydt[V.{species}] = ", + f"dydt[V.{species}] = 0 # ", + ) + for line_num, line in enumerate(lines): + if line.startswith(2 * self.indentation + "dydt = [0] * V.NUM\n"): + # Write right-hand side of the differential equation: dydt + lines[line_num + 1] = ( + 2 * self.indentation + + f"\n{2 * self.indentation}".join(self.differential_equations) + + "\n\n" + ) + elif line.startswith(self.indentation + "return x"): + if not self.param_info: + # Set all parameter values = 1.0 if param_info is not given + lines[line_num - 1] = ( + f"{self.indentation + 'x[C.'}" + + f"] = 1.0\n{self.indentation + 'x[C.'}".join(self.parameters) + + "] = 1.0\n\n" + ) + else: + # Set parameter values using param_info + lines[line_num - 1] = ( + "{}".format(self.indentation) + + "\n{}".format(self.indentation).join(self.param_info) + + "\n\n" + ) + elif line.startswith(self.indentation + "return y0"): + if self.init_info: + # Set initial values using init_info + lines[line_num - 1] = ( + "{}".format(self.indentation) + + "\n{}".format(self.indentation).join(self.init_info) + + "\n\n" + ) + with open( + os.path.join( + f"{self.name}", + "ode.py", + ), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + with open( + os.path.join(os.path.dirname(__file__), "template", "reaction_network.py"), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + for line_num, line in enumerate(lines): + if line.startswith(2 * self.indentation + "v = {}\n"): + # Write flux vector: v + lines[line_num + 1] = ( + 2 * self.indentation + + f"\n{2 * self.indentation}".join(self.reactions) + + "\n\n" + ) + with open( + os.path.join(f"{self.name}", "reaction_network.py"), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + else: + lines = jl.ODE.splitlines() + for line_num, line in enumerate(lines): + if line.startswith(self.indentation + "v = Dict{Int64,Float64}()"): + # Write flux vector: v + lines[line_num + 1] = ( + self.indentation + f"\n{self.indentation}".join(self.reactions) + "\n" + ) + lines[line_num + 1] = ( + lines[line_num + 1] + .replace("x[C.", "p[C.") + .replace("y[V.", "u[V.") + .replace("**", "^") + ) + # Write right-hand side of the differential equation: dydt + lines[line_num + 5] = ( + self.indentation + + f"\n{self.indentation}".join(self.differential_equations) + + "\n" + ).replace("dydt[V.", "du[V.") + elif line.startswith(self.indentation + "return p"): + if not self.param_info: + # Set all parameter values = 1.0 if param_info is not given + lines[line_num - 1] = ( + f"{self.indentation + 'p[C.'}" + + f"] = 1.0\n{self.indentation + 'p[C.'}".join(self.parameters) + + "] = 1.0\n\n" + ) + else: + # Set parameter values using param_info + lines[line_num - 1] = ( + "{}".format(self.indentation) + + "\n{}".format(self.indentation).join(self.param_info) + + "\n\n" + ).replace("x[C.", "p[C.") + elif line.startswith(self.indentation + "return u0"): + if self.init_info: + # Set initial values using init_info + lines[line_num - 1] = ( + "{}".format(self.indentation) + + "\n{}".format(self.indentation).join(self.init_info) + + "\n\n" + ).replace("y0[V.", "u0[V.") + with open( + os.path.join( + f"{self.name}_jl", + "ode.jl", + ), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + + def _update_search_param(self) -> None: + """ + Update search_param.py + """ + if self.lang == "python": + with open( + os.path.join(os.path.dirname(__file__), "template", "search_param.py"), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + for line_num, line in enumerate(lines): + if "self.idx_params = []" in line: + # Write all parameters in idx_params (except for param_excluded) + lines[line_num] = f"{2 * self.indentation}self.idx_params = [\n" + lines[line_num] += ( + "{}".format(3 * self.indentation + "C.") + + ",\n{}".format(3 * self.indentation + "C.").join(self.parameters) + + ",\n" + ) + lines[line_num] += f"{2 * self.indentation}]\n" + for param_name in self.param_excluded: + # Comment out parameters in param_excluded + lines[line_num] = lines[line_num].replace( + f"C.{param_name},", f"# C.{param_name}," + ) + elif "# parameter constraints" in line: + # Add parameter constraints + lines[line_num + 1] = ( + "{}".format(2 * self.indentation) + + "\n{}".format(2 * self.indentation).join(self.param_constraints) + + "\n\n" + ) + with open( + os.path.join(f"{self.name}", "search_param.py"), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + else: + lines = jl.SEARCH_PARAM.splitlines() + for line_num, line in enumerate(lines): + if "search_idx_params::Vector{Int} = []" in line: + # Write all parameters in idx_params (except for param_excluded) + lines[line_num] = ( + f"{self.indentation}search_idx_params" + "::Vector{Int} = [\n" + ) + lines[line_num] += ( + "{}".format(2 * self.indentation + "C.") + + ",\n{}".format(2 * self.indentation + "C.").join(self.parameters) + + ",\n" + ) + lines[line_num] += f"{self.indentation}]\n" + for param_name in self.param_excluded: + # Comment out parameters in param_excluded + lines[line_num] = lines[line_num].replace( + f"C.{param_name},", f"# C.{param_name}," + ) + elif "# parameter constraints" in line: + # Add parameter constraints + lines[line_num + 1] = ( + "{}".format(self.indentation) + + "\n{}".format(self.indentation).join(self.param_constraints) + + "\n" + ).replace("x[C.", "p[C.") + with open( + os.path.join(f"{self.name}_jl", "search_param.jl"), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + + def _convert_names( + self, + line: str, + *, + p: Optional[List[str]] = None, + u: Optional[List[str]] = None, + init: Optional[List[str]] = None, + ) -> str: + """ + Replace + - p[xxx] with x[C.xxx] + - u[xxx] with sol.y[V.xxx] + - init[xxx] with y0[V.xxx] + + Parameters + ---------- + p : list of strings, optional + Parameters. + + u : list of strings, optional + Species. + + init : list of strings, optional + Initial conditions. + + Returns + ------- + line : string + Each line. + + """ + if p is None: + p = [] + if u is None: + u = [] + if init is None: + init = [] + for p_name in p: + if not p_name.strip() in self.parameters: + raise NameError(f"{p_name.strip()} is not defined in model parameters.") + else: + line = ( + line.replace(f"p[{p_name.strip()}]", f"x[C.{p_name.strip()}]") + if self.lang == "python" + else line.replace(f"p[{p_name.strip()}]", f"p[C.{p_name.strip()}]") + ) + for s_name in u: + if not s_name.strip() in self.species: + raise NameError(f"{s_name.strip()} is not defined in model species.") + else: + line = ( + line.replace(f"u[{s_name.strip()}]", f"sol.y[V.{s_name.strip()}]") + if self.lang == "python" + else line.replace(f"u[{s_name.strip()}]", f"sol.u[j][V.{s_name.strip()}]") + ) + for s_name in init: + if not s_name.strip() in self.species: + raise NameError(f"{s_name.strip()} is not defined in model species.") + else: + line = ( + line.replace(f"init[{s_name.strip()}]", f"y0[V.{s_name.strip()}]") + if self.lang == "python" + else line.replace(f"init[{s_name.strip()}]", f"u0[V.{s_name.strip()}]") + ) + return line + + def _update_observable(self) -> None: + """ + Update observable.py + """ + if self.lang == "python": + with open( + os.path.join(os.path.dirname(__file__), "template", "observable.py"), + encoding="utf-8", + mode="r", + ) as f: + lines = f.readlines() + if not self.sim_conditions: + # When sim condition is not given, add 'sim condition control: pass' + self.sim_conditions = [["control", "pass"]] + for line_num, line in enumerate(lines): + if line.startswith(f"{2 * self.indentation}self.obs_names: list = []"): + # Write observables + lines[line_num] = f"{2 * self.indentation}self.obs_names: list = [\n" + lines[line_num] += ( + f"{3 * self.indentation}'" + + f"',\n{3 * self.indentation}'".join( + [desc[0].strip() for desc in self.obs_desc] + ) + + "',\n" + ) + lines[line_num] += f"{2 * self.indentation}]\n" + elif f"self.t: range = range(101)" in line: + # Write interval of integration + if self.sim_tspan: + lines[line_num] = "{}self.t: range = range({}, {}+1)\n".format( + 2 * self.indentation, + self.sim_tspan[0].strip(), + self.sim_tspan[1].strip(), + ) + elif line.startswith(f"{2 * self.indentation}self.conditions: list = []"): + # Write sim.condition + lines[line_num] = f"{2 * self.indentation}self.conditions: list = [\n" + lines[line_num] += ( + f"{3 * self.indentation}'" + + f"',\n{3 * self.indentation}'".join( + [condition[0].strip() for condition in self.sim_conditions] + ) + + "',\n" + ) + lines[line_num] += f"{2 * self.indentation}]\n" + elif "# unperturbed steady state" in line: + if self.sim_unperturbed: + lines[line_num + 1] = ( + 2 * self.indentation + + f"\n{2 * self.indentation}".join( + c.strip() for c in self.sim_unperturbed.split(sep=";") + ) + + "\n" + ) + lines[line_num + 1] += ( + 2 * self.indentation + + "y0 = get_steady_state(self.diffeq, y0, tuple(x))\n" + + 2 * self.indentation + + f"if not y0:\n{3 * self.indentation}return False\n" + ) + # pa: parameters + # init: initial conditions + lines[line_num + 1] = self._convert_names( + line=lines[line_num + 1], + p=re.findall(r"p\[(.*?)\]", self.sim_unperturbed), + init=re.findall(r"init\[(.*?)\]", self.sim_unperturbed), + ) + elif "for i, condition in enumerate(self.conditions):" in line: + for i, condition in enumerate(self.sim_conditions): + # Use ';' for adding description of each condition + if i == 0: + lines[line_num + 1] = ( + 3 * self.indentation + + f'if condition == \'{condition[0].strip(" ")}\':\n' + + 4 * self.indentation + + f"\n{4 * self.indentation}".join( + c.strip(" ") for c in condition[1].split(sep=";") + ) + + ("\n\n" if len(self.sim_conditions) == 1 else "") + ) + else: + lines[line_num + 1] += ( + 3 * self.indentation + + f'elif condition == \'{condition[0].strip(" ")}\':\n' + + 4 * self.indentation + + f"\n{4 * self.indentation}".join( + c.strip(" ") for c in condition[1].split(sep=";") + ) + + ("\n\n" if i == len(self.sim_conditions) - 1 else "") + ) + # pa: parameters + # init: initial conditions + lines[line_num + 1] = self._convert_names( + line=lines[line_num + 1], + p=re.findall(r"p\[(.*?)\]", condition[1]), + init=re.findall(r"init\[(.*?)\]", condition[1]), + ) + elif "sol is None:" in line: + lines[line_num + 3] = "" # initialization (default: pass) + for desc in self.obs_desc: + lines[line_num + 3] += ( + "{}".format(4 * self.indentation) + + "self.simulations" + + f"[self.obs_names.index('{desc[0].strip()}'), i] = (\n" + + f"{5 * self.indentation}" + + desc[1].strip(" ").strip() + ) + # p: parameters + # u: species + lines[line_num + 3] = self._convert_names( + line=lines[line_num + 3], + p=re.findall(r"p\[(.*?)\]", desc[1]), + u=re.findall(r"u\[(.*?)\]", desc[1]), + ) + lines[line_num + 3] += "\n{}".format(4 * self.indentation + ")\n") + with open( + os.path.join(f"{self.name}", "observable.py"), + encoding="utf-8", + mode="w", + ) as f: + f.writelines(lines) + else: + self._split_observables_julia() + + def _split_observables_julia(self) -> None: + """ + Create observable.jl, simulation.jl and experimental_data.jl + """ + # observable.jl + lines = jl.OBSERVABLE.splitlines() + for line_num, line in enumerate(lines): + if "const observables = []" in line: + lines[line_num] = "const observables = [\n" + lines[line_num + 1] = ( + f'{self.indentation}"' + + f'",\n{self.indentation}"'.join([desc[0].strip() for desc in self.obs_desc]) + + '",\n' + ) + lines[line_num + 1] += "]\n" + with open( + os.path.join(f"{self.name}_jl", "observable.jl"), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + # simulation.jl + lines = jl.SIMULATION.splitlines() + for line_num, line in enumerate(lines): + if line.startswith("const t = collect(0.0:dt:100.0)"): + # Write interval of integration + if self.sim_tspan: + lines[line_num] = "const t = collect({}:dt:{})".format( + self.sim_tspan[0].strip(), + self.sim_tspan[1].strip(), + ) + elif line.startswith("const conditions = []"): + # Write sim.condition + lines[line_num] = "const conditions = [" + lines[line_num + 1] = ( + f'{self.indentation}"' + + f'",\n{self.indentation}"'.join( + [condition[0].strip() for condition in self.sim_conditions] + ) + + '",\n' + ) + lines[line_num + 1] += "]\n" + elif "# unperturbed steady state" in line: + if self.sim_unperturbed: + lines[line_num + 1] = ( + self.indentation + + f"\n{self.indentation}".join( + c.strip() for c in self.sim_unperturbed.split(sep=";") + ) + + "\n" + ) + lines[line_num + 1] += ( + f"{self.indentation}u0 = get_steady_state(diffeq!, u0, p)\n" + + f"{self.indentation}if isempty(u0)\n" + + f"{2 * self.indentation}return false\n" + + f"{self.indentation}end\n" + ) + # p: parameters + # init: initial conditions + lines[line_num + 1] = self._convert_names( + line=lines[line_num + 1], + p=re.findall(r"p\[(.*?)\]", self.sim_unperturbed), + init=re.findall(r"init\[(.*?)\]", self.sim_unperturbed), + ) + elif "for (i, condition) in enumerate(conditions)" in line: + for i, condition in enumerate(self.sim_conditions): + # Use ';' for adding description of each condition + if i == 0: + lines[line_num + 1] = ( + 3 * self.indentation + + f'if condition == "{condition[0].strip(" ")}"\n' + + 4 * self.indentation + + f"\n{4 * self.indentation}".join( + c.strip(" ") for c in condition[1].split(sep=";") + ) + + ("\n\n" if len(self.sim_conditions) == 1 else "") + ) + else: + lines[line_num + 1] += ( + 3 * self.indentation + + f'elseif condition == "{condition[0].strip(" ")}"\n' + + 4 * self.indentation + + f"\n{4 * self.indentation}".join( + c.strip(" ") for c in condition[1].split(sep=";") + ) + + f"\n{3 * self.indentation}end\n" + ) + # p: parameters + # init: initial conditions + lines[line_num + 1] = self._convert_names( + line=lines[line_num + 1], + p=re.findall(r"p\[(.*?)\]", condition[1]), + init=re.findall(r"init\[(.*?)\]", condition[1]), + ) + elif "if sol === nothing" in line: + lines[line_num + 4] = "" # initialization (default: # line_num + 4) + for desc in self.obs_desc: + lines[line_num + 4] += ( + "{}".format(4 * self.indentation) + + "simulations" + + f'[observables_index("{desc[0].strip()}"), i, j] = (\n' + + f"{5 * self.indentation}" + + desc[1].strip(" ").strip() + ) + # p: parameters + # u: species + lines[line_num + 4] = self._convert_names( + line=lines[line_num + 4], + p=re.findall(r"p\[(.*?)\]", desc[1]), + u=re.findall(r"u\[(.*?)\]", desc[1]), + ) + lines[line_num + 4] += "\n{}".format(4 * self.indentation + ")\n") + with open( + os.path.join(f"{self.name}_jl", "simulation.jl"), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + # experimental_data.jl + lines = jl.EXPERIMENTAL_DATA.splitlines() + with open( + os.path.join(f"{self.name}_jl", "experimental_data.jl"), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + + def convert( + self, + *, + show_restrictions: bool = False, + overwrite: bool = False, + ) -> None: + """ + Convert text to a biomass-formatted model. + + Parameters + ---------- + show_restrictions : bool (defauld: :obj:`False`) + Whether to display reaction indices in which thermodynamic restrictions should be + imposed. These detailed balance constraints require the product of the equilibrium + constants along a cycle to be equal to 1. + overwrite : bool (defauld: :obj:`False`) + If :obj:`True`, the model folder will be overwritten. + + Examples + -------- + >>> from pasmopy import Text2Model + >>> Text2Model("Kholodenko1999.txt").convert() + + """ + if overwrite and os.path.isdir( + f"{self.name}" if self.lang == "python" else f"{self.name}_jl" + ): + shutil.rmtree(f"{self.name}" if self.lang == "python" else f"{self.name}_jl") + if self.lang == "python": + shutil.copytree( + os.path.join( + os.path.dirname(__file__), + "template", + ), + f"{self.name}", + ) + else: + os.makedirs(os.path.join(f"{self.name}_jl", "name2idx")) + + self.create_ode() + self.find_cyclic_reaction_routes() + self._update_parameters() + self._update_species() + self._update_diffeq() + self._update_search_param() + self._update_observable() + if self.lang == "julia": + # Create fitness.jl + lines = jl.PROBLEM.splitlines() + with open( + os.path.join(f"{self.name}_jl", "problem.jl"), + encoding="utf-8", + mode="w", + ) as f: + f.write("\n".join(lines)) + print("Model information\n-----------------") + print(f"{len(self.reactions):d} reactions") + print(f"{len(self.species):d} species") + print(f"{len(self.parameters):d} parameters") + if show_restrictions: + if len(self.restrictions) == 0: + print("No cyclic reaction routes.") + else: + print("\nThermodynamic restrictions") + print("--------------------------") + for restriction in self.restrictions: + print("{" + ", ".join(restriction) + "}") + + def to_markdown(self, num_reactions: int, savedir: str = "markdown") -> None: + """ + Create markdown table describing differential equations. + + Parameters + ---------- + num_reactions : int + The number of rate equations in the model. + savedir : str (default: "markdown") + The directory name to save the output. + + Examples + -------- + >>> from pasmopy import Text2Model + >>> Text2Model("Kholodenko1999.txt").to_markdown(25) + + """ + os.makedirs(os.path.join(savedir, f"{self.name.split(os.sep)[-1]}"), exist_ok=True) + self.create_ode() + with open(self.input_txt, encoding="utf-8") as f: + lines = f.readlines() + for num, line in enumerate(lines): + if num_reactions <= num: + break + else: + rate_equation_formatted = f"{self.reactions[num]}".split("=")[1].strip() + for p_names in self.parameters: + if f" ** x[C.{p_names.strip()}]" in f"{self.reactions[num]}": + # Hill coefficients + rate_equation_formatted = rate_equation_formatted.replace( + f" ** x[C.{p_names.strip()}]", + f"{p_names.strip()}".replace( + f"{num + 1:d}", f"{num + 1:d}" + ), + ) + elif f"x[C.{p_names.strip()}]" in f"{self.reactions[num]}": + # Others + rate_equation_formatted = rate_equation_formatted.replace( + f"x[C.{p_names.strip()}]", + f"{p_names.strip()}".replace( + f"{num + 1:d}", f"{num + 1:d}" + ), + ) + reaction = line.split("|")[0].rstrip() + if reaction.startswith("@rxn "): + reaction = self._remove_prefix(reaction, "@rxn ").split(":")[0].strip() + lines[num] = ( + f"|{num + 1:d}|" + + reaction + + f"|{rate_equation_formatted.replace('*', '·').replace('y[V.', '[')}|" + + "\n" + ) + with open( + os.path.join(savedir, f"{self.name.split(os.sep)[-1]}", "rate.md"), + encoding="utf-8", + mode="w", + ) as f: + f.writelines( + [ + "|No.|Description|Rate|\n|---|-----------|----|\n", + *lines[:num_reactions], + ] + ) + # Differential equation + differential_equations_formatted = [ + eq.replace("*", "·") for eq in self.differential_equations + ] + for i, eq in enumerate(differential_equations_formatted): + for s_name in self.species: + if f"dydt[V.{s_name.strip()}]" in eq: + eq = eq.replace( + f"dydt[V.{s_name.strip()}]", + f"|{i + 1:d}|d[{s_name.strip()}]/dt", + ) + break + for num in range(num_reactions): + if f"v[{num + 1:d}]" in eq: + eq = eq.replace( + f"v[{num + 1:d}]", + f"_v_ {num + 1:d}", + ) + differential_equations_formatted[i] = eq + "|\n" + with open( + os.path.join(savedir, f"{self.name.split(os.sep)[-1]}", "diffeq.md"), + encoding="utf-8", + mode="w", + ) as f: + f.writelines( + [ + "|No.|Differential equations|\n|---|----------------------|\n", + *differential_equations_formatted, + ] + ) + + def register_word( + self, + terminology: Optional[ + Dict[ + Literal[ + "_bind_and_dissociate", + "dimerize", + "bind", + "dissociate", + "is_phosphorylated", + "is_dephosphorylated", + "phosphorylate", + "dephosphorylate", + "transcribe", + "is_translated", + "synthesize", + "is_synthesized", + "degrade", + "is_degraded", + "translocate", + ], + List[str], + ] + ] = None, + ) -> None: + """ + Register user-defined rule word. + + Parameters + ---------- + terminology : Dict[str, List[str]], optional + Rule to which users register a new rule word and list of user-defined rule words. + + Examples + -------- + >>> from pasmopy import Text2Model + >>> mm_kinetics = Text2Model("michaelis_menten.txt") + >>> mm_kinetics.register_word({"dissociate": ["releases"]}) + >>> mm_kinetics.convert() + + """ + if terminology is None: + terminology = {} + for rxn_rule in terminology.keys(): + assert isinstance(rxn_rule, str) + if rxn_rule not in self.rule_words.keys(): + raise ValueError( + f"{rxn_rule} is not defined in reaction_rules.\n" + f"Choose a reaction rule from {', '.join(map(str, self.rule_words.keys()))}" + ) + for my_word in terminology[rxn_rule]: + assert isinstance(my_word, str) + for rule, words in self.rule_words.items(): + for registered_word in words: + if " " + my_word in registered_word and registered_word in " " + my_word: + raise NameError( + f"Cannot supply '{my_word}' to '{rxn_rule}'. " + f"Currently, it is used in the rule: {rule}" + ) + self.rule_words[rxn_rule].append(" " + my_word) diff --git a/biomass/construction/thermodynamic_restrictions.py b/biomass/construction/thermodynamic_restrictions.py new file mode 100644 index 00000000..3ce5363c --- /dev/null +++ b/biomass/construction/thermodynamic_restrictions.py @@ -0,0 +1,205 @@ +from collections import Counter +from dataclasses import dataclass, field +from typing import List, NamedTuple, Tuple + + +class DuplicateError(Exception): + pass + + +class ComplexFormation(NamedTuple): + rxn_no: int + components: set + complex: str + is_binding: bool + + +@dataclass +class ThermodynamicRestrictions(object): + """ + Thermodynamic restrictions along cyclic pathways in a kinetic scheme. + + Notes + ----- + If a kinetic scheme includes "true" cycles, in which the initial and final states are + identical, the equilibrium constants of the reactions along any cycle satisfy so-called + "detailed balance" relationships. These detailed balance relations require the product + of the equilibrium constants along a cycle to be equal to 1, since at equilibrium the net flux + through any cycle vanishes. + """ + + _rxn_indices: dict = field( + default_factory=dict, + init=False, + ) + _tree: dict = field( + default_factory=dict, + init=False, + ) + restrictions: list = field( + default_factory=list, + init=False, + ) + complex_formations: List[ComplexFormation] = field( + default_factory=list, + init=False, + ) + + def _create_tree( + self, + tree: dict, + parent_pattern: ComplexFormation, + ) -> dict: + """ + Create tree. + """ + for component in parent_pattern.components: + for another_pattern in self.complex_formations: + if ( + component == another_pattern.complex + and parent_pattern.rxn_no != another_pattern.rxn_no + ): + for k in another_pattern.components: + if component not in tree.keys(): + tree.setdefault(component, {}) + tree[component].setdefault(k, {f"{another_pattern.rxn_no}": {}}) + tree[component][k] = self._create_tree(tree[component][k], another_pattern) + return tree + + def _add_to_rxn_indices1(self, complex_name: str, tree: dict) -> None: + for component in tree.keys(): + if not tree[component]: + self._rxn_indices[complex_name].append(component) + else: + self._add_to_rxn_indices1(complex_name, tree[component]) + + def _add_to_rxn_indices2(self, complex_name: str, monomer: str, pair: str) -> None: + for another_complex in self._tree[complex_name].keys(): + if ( + monomer in self._tree[complex_name][another_complex].keys() + and len(self._tree[complex_name][another_complex][monomer]) > 1 + ): + for name in self._tree[complex_name][another_complex][monomer].keys(): + if name in self._tree[complex_name][pair].keys(): + self._append_reaction_number( + complex_name, + monomer, + pair, + name, + another_complex, + ) + + def _add_to_rxn_indices3(self, complex_name: str, monomer: str, pair: str) -> None: + for another_complex in self._tree[complex_name].keys(): + if ( + monomer in self._tree[complex_name][another_complex].keys() + and len(self._tree[complex_name][another_complex][monomer]) == 1 + ): + for name in self._tree[complex_name][another_complex].keys(): + if not name.isdecimal() and name != monomer: + for components in self._tree[complex_name][another_complex].keys(): + if name == components: + for k in self._tree[complex_name][pair][name].keys(): + if not k.isdecimal() and k in self._tree[complex_name].keys(): + self._append_reaction_number( + complex_name, + monomer, + pair, + name, + another_complex, + ) + + def _append_reaction_number( + self, + complex_name: str, + monomer: str, + pair: str, + name: str, + another_complex: str, + ) -> None: + for val in self._tree[complex_name][another_complex].keys(): + if val.isdecimal(): + self._rxn_indices[monomer].append(val) + for val in self._tree[complex_name][another_complex][monomer].keys(): + if val.isdecimal(): + self._rxn_indices[monomer].append(val) + for val in self._tree[complex_name][pair][name].keys(): + if val.isdecimal(): + self._rxn_indices[monomer].append(val) + + def _get_complex_patterns(self) -> List[Tuple[ComplexFormation, ComplexFormation]]: + complex_patterns = [] + for i, pattern_a in enumerate(self.complex_formations): + for j in range(i + 1, len(self.complex_formations)): + pattern_b = self.complex_formations[j] + if ( + pattern_a.components == pattern_b.components + and pattern_a.complex == pattern_b.complex + ): + raise DuplicateError( + "Duplicate binding-dissociation events are detected " + f"in lines {pattern_a.rxn_no:d} and {pattern_b.rxn_no:d}." + ) + elif ( + pattern_a.is_binding != pattern_b.is_binding + and pattern_a.complex == pattern_b.complex + ): + complex_patterns.append((pattern_a, pattern_b)) + return complex_patterns + + def _impose_restrictions(self, complex_name: str) -> None: + count = Counter(self._rxn_indices[complex_name]) + if all(n == 2 for n in count.values()): + self.restrictions.append(list(set(self._rxn_indices[complex_name]))) + else: + _monomers = {} + for name in self._tree[complex_name].keys(): + if len(self._tree[complex_name][name].values()) == 1: + _monomers[name] = list(self._tree[complex_name][name].keys())[0] + for monomer, ridx in _monomers.items(): + self._rxn_indices[monomer] = [ridx] + pair = None + for another_complex in self._tree[complex_name].keys(): + if ( + another_complex != monomer + and ridx in self._tree[complex_name][another_complex] + ): + pair = another_complex + if pair is not None: + self._add_to_rxn_indices2(complex_name, monomer, pair) + _reactions = list(set(self._rxn_indices[monomer])) + if len(_reactions) > 2 and _reactions not in self.restrictions: + self.restrictions.append(_reactions) + # initialize monomers + for monomer, ridx in _monomers.items(): + self._rxn_indices[monomer] = [ridx] + pair = None + for another_complex in self._tree[complex_name].keys(): + if ( + another_complex != monomer + and ridx in self._tree[complex_name][another_complex] + ): + pair = another_complex + if pair is not None: + self._add_to_rxn_indices3(complex_name, monomer, pair) + _reactions = list(set(self._rxn_indices[monomer])) + if len(_reactions) > 2 and _reactions not in self.restrictions: + self.restrictions.append(_reactions) + + def find_cyclic_reaction_routes(self) -> None: + """ + Find cyclic pathways in a reaction network. + """ + complex_patterns = self._get_complex_patterns() + _tree = {} + for patterns in complex_patterns: + for pattern in patterns: + _tree.setdefault(pattern.complex, {}) + for component in pattern.components: + _tree[pattern.complex].setdefault(component, {f"{pattern.rxn_no}": {}}) + _tree[pattern.complex] = self._create_tree(_tree[pattern.complex], pattern) + self._tree = _tree + for complex_name in self._tree.keys(): + self._rxn_indices.setdefault(complex_name, []) + self._add_to_rxn_indices1(complex_name, self._tree[complex_name]) + self._impose_restrictions(complex_name) diff --git a/docs/api/construction.rst b/docs/api/construction.rst new file mode 100644 index 00000000..72b675d7 --- /dev/null +++ b/docs/api/construction.rst @@ -0,0 +1,10 @@ +Construction of executable models from text (:py:mod:`biomass.construction`) +============================================================================ + +.. toctree:: + :maxdepth: 3 + + thermodynamic_restrictions + reaction_rules + text2model + template \ No newline at end of file diff --git a/docs/api/index.rst b/docs/api/index.rst index 4c4d18d9..d4a1d11c 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -6,6 +6,7 @@ API graph model_object + construction core result optimizer diff --git a/docs/api/reaction_rules.rst b/docs/api/reaction_rules.rst new file mode 100644 index 00000000..8b15e75a --- /dev/null +++ b/docs/api/reaction_rules.rst @@ -0,0 +1,7 @@ +Available reaction rules (:py:mod:`biomass.construction.reaction_rules`) +------------------------------------------------------------------------ + +.. autoclass:: biomass.construction.reaction_rules.ReactionRules + :members: _bind_and_dissociate, dimerize, bind, is_dissociated, is_phosphorylated, is_dephosphorylated, phosphorylate, dephosphorylate, transcribe, synthesize, is_synthesized, degrade, is_degraded, translocate, user_defined + :member-order: bysource + diff --git a/docs/api/template.rst b/docs/api/template.rst new file mode 100644 index 00000000..e0ed22c4 --- /dev/null +++ b/docs/api/template.rst @@ -0,0 +1,20 @@ +Template for BioMASS model construction (:py:mod:`biomass.construction.template`) +--------------------------------------------------------------------------------- + +.. automodule:: biomass.construction.template.ode + :members: + +.. automodule:: biomass.construction.template.observable + :members: + +.. automodule:: biomass.construction.template.search_param + :members: + +.. automodule:: biomass.construction.template.problem + :members: + +.. automodule:: biomass.construction.template.viz + :members: + +.. automodule:: biomass.construction.template.reaction_network + :members: \ No newline at end of file diff --git a/docs/api/text2model.rst b/docs/api/text2model.rst new file mode 100644 index 00000000..fa49f4a4 --- /dev/null +++ b/docs/api/text2model.rst @@ -0,0 +1,6 @@ +Text-to-model conversion (:py:mod:`biomass.construction.text2model`) +-------------------------------------------------------------------- + +.. autoclass:: biomass.construction.text2model.Text2Model + :members: convert, to_markdown, register_word + diff --git a/docs/api/thermodynamic_restrictions.rst b/docs/api/thermodynamic_restrictions.rst new file mode 100644 index 00000000..3ba191ef --- /dev/null +++ b/docs/api/thermodynamic_restrictions.rst @@ -0,0 +1,6 @@ +Thermodynamic restrictions (:py:mod:`biomass.construction.thermodynamic_restrictions`) +-------------------------------------------------------------------------------------- + +.. autoclass:: biomass.construction.thermodynamic_restrictions.ThermodynamicRestrictions + :members: + diff --git a/tests/test_text2model/C.py b/tests/test_text2model/C.py new file mode 100755 index 00000000..a4b9a691 --- /dev/null +++ b/tests/test_text2model/C.py @@ -0,0 +1,235 @@ +INDENT = " " * 4 + + +PPMEK = """\ + + @staticmethod + def _get_ppMEK_slope(t, ligand): + assert ligand in ['EGF', 'HRG'] + timepoints = [0, 300, 600, 900, 1800, 2700, 3600, 5400] + ppMEK_data = { + 'EGF': [0.000, 0.773, 0.439, 0.252, 0.130, 0.087, 0.080, 0.066], + 'HRG': [0.000, 0.865, 1.000, 0.837, 0.884, 0.920, 0.875, 0.789], + } + assert len(timepoints) == len(ppMEK_data[ligand]) + slope = [ + (ppMEK_data[ligand][i + 1] - activity) / (timepoints[i + 1] - timepoint) + for i, (timepoint, activity) in enumerate(zip(timepoints, ppMEK_data[ligand])) + if i + 1 < len(timepoints) + ] + for i, timepoint in enumerate(timepoints): + if timepoint <= t <= timepoints[i + 1]: + return slope[i] + assert False + +""" + + +CONDITIONS = """\ + + if x[C.Ligand] == 1: # EGF=10nM + dydt[V.ppMEKc] = self._get_ppMEK_slope(t, 'EGF') + elif x[C.Ligand] == 2: # HRG=10nM + dydt[V.ppMEKc] = self._get_ppMEK_slope(t, 'HRG') + else: # Default: No ligand input + dydt[V.ppMEKc] = 0.0 + +""" + + +BOUNDS = """\ + + search_rgn[:, C.V1] = [7.33e-2, 6.60e-01] + search_rgn[:, C.K1] = [1.83e2, 8.50e2] + search_rgn[:, C.V5] = [6.48e-3, 7.20e1] + search_rgn[:, C.K5] = [6.00e-1, 1.60e04] + search_rgn[:, C.V10] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K10] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.n10] = [1.00, 4.00] + search_rgn[:, C.kf11] = [8.30e-13, 1.44e-2] + search_rgn[:, C.kf12] = [8.00e-8, 5.17e-2] + search_rgn[:, C.kf13] = [1.38e-7, 4.84e-1] + search_rgn[:, C.V14] = [4.77e-3, 4.77e1] + search_rgn[:, C.K14] = [2.00e2, 2.00e6] + search_rgn[:, C.V15] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K15] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kf18] = [2.20e-4, 5.50e-1] + search_rgn[:, C.kr18] = [2.60e-4, 6.50e-1] + search_rgn[:, C.V20] = [4.77e-3, 4.77e1] + search_rgn[:, C.K20] = [2.00e2, 2.00e6] + search_rgn[:, C.V21] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K21] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.V24] = [4.77e-2, 4.77e0] + search_rgn[:, C.K24] = [2.00e3, 2.00e5] + search_rgn[:, C.V25] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K25] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kf26] = [2.20e-4, 5.50e-1] + search_rgn[:, C.kr26] = [2.60e-4, 6.50e-1] + search_rgn[:, C.V27] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K27] = [1.00e2, 1.00e4] + search_rgn[:, C.V28] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K28] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.V29] = [4.77e-2, 4.77e0] + search_rgn[:, C.K29] = [2.93e3, 2.93e5] + search_rgn[:, C.V30] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K30] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.V31] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K31] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.n31] = [1.00, 4.00] + search_rgn[:, C.kf32] = [8.30e-13, 1.44e-2] + search_rgn[:, C.kf33] = [8.00e-8, 5.17e-2] + search_rgn[:, C.kf34] = [1.38e-7, 4.84e-1] + search_rgn[:, C.V35] = [4.77e-3, 4.77e1] + search_rgn[:, C.K35] = [2.00e2, 2.00e6] + search_rgn[:, C.V36] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K36] = [1.00e2, 1.00e4] + search_rgn[:, C.V37] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K37] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kf40] = [2.20e-4, 5.50e-1] + search_rgn[:, C.kr40] = [2.60e-4, 6.50e-1] + search_rgn[:, C.V42] = [4.77e-3, 4.77e1] + search_rgn[:, C.K42] = [2.00e2, 2.00e6] + search_rgn[:, C.V43] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K43] = [1.00e2, 1.00e4] + search_rgn[:, C.V44] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K44] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kf47] = [1.45e-4, 1.45e0] + search_rgn[:, C.kr47] = [6.00e-3, 6.00e1] + search_rgn[:, C.kf48] = [2.70e-3, 2.70e1] + search_rgn[:, C.kf49] = [5.00e-5, 5.00e-1] + search_rgn[:, C.kr49] = [5.00e-3, 5.00e1] + search_rgn[:, C.kf50] = [3.00e-3, 3.00e1] + search_rgn[:, C.kf51] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kr51] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.V57] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.K57] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.n57] = [1.00, 4.00] + search_rgn[:, C.kf58] = [8.30e-13, 1.44e-2] + search_rgn[:, C.kf59] = [8.00e-8, 5.17e-2] + search_rgn[:, C.kf60] = [1.38e-7, 4.84e-1] + search_rgn[:, C.kf61] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.kf62] = [2.20e-4, 5.50e-1] + search_rgn[:, C.kr62] = [2.60e-4, 6.50e-1] + search_rgn[:, C.kf63] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.KF31] = [np.exp(-10), np.exp(10)] + search_rgn[:, C.nF31] = [1.00, 4.00] + search_rgn[:, C.a] = [1.00e2, 5.00e2] + +""" + +NORMALIZATION = """\ + for observable in self.obs_names: + self.normalization[observable] = {"timepoint": None, "condition": []} +""" + +EXPERIMENTAL_DATA = """\ + def set_data(self): + self.experiments[self.obs_names.index("Phosphorylated_MEKc")] = { + "EGF": [0.000, 0.773, 0.439, 0.252, 0.130, 0.087, 0.080, 0.066], + "HRG": [0.000, 0.865, 1.000, 0.837, 0.884, 0.920, 0.875, 0.789], + } + self.error_bars[self.obs_names.index("Phosphorylated_MEKc")] = { + "EGF": [ + sd / np.sqrt(3) for sd in [0.000, 0.030, 0.048, 0.009, 0.009, 0.017, 0.012, 0.008] + ], + "HRG": [ + sd / np.sqrt(3) for sd in [0.000, 0.041, 0.000, 0.051, 0.058, 0.097, 0.157, 0.136] + ], + } + + self.experiments[self.obs_names.index("Phosphorylated_ERKc")] = { + "EGF": [0.000, 0.867, 0.799, 0.494, 0.313, 0.266, 0.200, 0.194], + "HRG": [0.000, 0.848, 1.000, 0.971, 0.950, 0.812, 0.747, 0.595], + } + self.error_bars[self.obs_names.index("Phosphorylated_ERKc")] = { + "EGF": [ + sd / np.sqrt(3) for sd in [0.000, 0.137, 0.188, 0.126, 0.096, 0.087, 0.056, 0.012] + ], + "HRG": [ + sd / np.sqrt(3) for sd in [0.000, 0.120, 0.000, 0.037, 0.088, 0.019, 0.093, 0.075] + ], + } + + self.experiments[self.obs_names.index("Phosphorylated_RSKw")] = { + "EGF": [0, 0.814, 0.812, 0.450, 0.151, 0.059, 0.038, 0.030], + "HRG": [0, 0.953, 1.000, 0.844, 0.935, 0.868, 0.779, 0.558], + } + self.error_bars[self.obs_names.index("Phosphorylated_RSKw")] = { + "EGF": [ + sd / np.sqrt(3) for sd in [0, 0.064, 0.194, 0.030, 0.027, 0.031, 0.043, 0.051] + ], + "HRG": [ + sd / np.sqrt(3) for sd in [0, 0.230, 0.118, 0.058, 0.041, 0.076, 0.090, 0.077] + ], + } + + self.experiments[self.obs_names.index("Phosphorylated_cFos")] = { + "EGF": [0, 0.060, 0.109, 0.083, 0.068, 0.049, 0.027, 0.017], + "HRG": [0, 0.145, 0.177, 0.158, 0.598, 1.000, 0.852, 0.431], + } + self.error_bars[self.obs_names.index("Phosphorylated_cFos")] = { + "EGF": [ + sd / np.sqrt(3) for sd in [0, 0.003, 0.021, 0.013, 0.016, 0.007, 0.003, 0.002] + ], + "HRG": [ + sd / np.sqrt(3) for sd in [0, 0.010, 0.013, 0.001, 0.014, 0.000, 0.077, 0.047] + ], + } + + # ---------------------------------------------------------------------- + + self.experiments[self.obs_names.index("Phosphorylated_CREBw")] = { + "EGF": [0, 0.446, 0.030, 0.000, 0.000], + "HRG": [0, 1.000, 0.668, 0.460, 0.340], + } + self.error_bars[self.obs_names.index("Phosphorylated_CREBw")] = { + "EGF": [sd / np.sqrt(3) for sd in [0, 0.0, 0.0, 0.0, 0.0]], + "HRG": [sd / np.sqrt(3) for sd in [0, 0.0, 0.0, 0.0, 0.0]], + } + # ---------------------------------------------------------------------- + + self.experiments[self.obs_names.index("cfos_mRNA")] = { + "EGF": [0, 0.181, 0.476, 0.518, 0.174, 0.026, 0.000], + "HRG": [0, 0.353, 0.861, 1.000, 0.637, 0.300, 0.059], + } + self.error_bars[self.obs_names.index("cfos_mRNA")] = { + "EGF": [sd / np.sqrt(3) for sd in [0.017, 0.004, 0.044, 0.004, 0.023, 0.007, 0.008]], + "HRG": [sd / np.sqrt(3) for sd in [0.017, 0.006, 0.065, 0.044, 0.087, 0.023, 0.001]], + } + # ---------------------------------------------------------------------- + + self.experiments[self.obs_names.index("cFos_Protein")] = { + "EGF": [0, 0.078, 0.216, 0.240, 0.320, 0.235], + "HRG": [0, 0.089, 0.552, 0.861, 1.000, 0.698], + } + self.error_bars[self.obs_names.index("cFos_Protein")] = { + "EGF": [sd / np.sqrt(3) for sd in [0, 0.036, 0.028, 0.056, 0.071, 0.048]], + "HRG": [sd / np.sqrt(3) for sd in [0, 0.021, 0.042, 0.063, 0.000, 0.047]], + } + + self.experiments[self.obs_names.index("dusp_mRNA")] = { + "EGF": [0.000, 0.177, 0.331, 0.214, 0.177, 0.231], + "HRG": [0.000, 0.221, 0.750, 1.000, 0.960, 0.934], + } + self.error_bars[self.obs_names.index("dusp_mRNA")] = { + "EGF": [sd / np.sqrt(3) for sd in [0.033, 0.060, 0.061, 0.032, 0.068, 0.050]], + "HRG": [sd / np.sqrt(3) for sd in [0.027, 0.059, 0.094, 0.124, 0.113, 0.108]], + } + + @staticmethod + def get_timepoint(obs_name) -> List[int]: + if obs_name in [ + "Phosphorylated_MEKc", + "Phosphorylated_ERKc", + "Phosphorylated_RSKw", + "Phosphorylated_cFos", + ]: + return [0, 300, 600, 900, 1800, 2700, 3600, 5400] # (Unit: sec.) + elif obs_name == "Phosphorylated_CREBw": + return [0, 600, 1800, 3600, 5400] + elif obs_name == "cfos_mRNA": + return [0, 600, 1200, 1800, 2700, 3600, 5400] + elif obs_name in ["cFos_Protein", "dusp_mRNA"]: + return [0, 900, 1800, 2700, 3600, 5400] + assert False +""" diff --git a/tests/test_text2model/__init__.py b/tests/test_text2model/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_text2model/simulations_original_BN.npy b/tests/test_text2model/simulations_original_BN.npy new file mode 100644 index 0000000000000000000000000000000000000000..bfdefa25b0a0df3e10542376b644489e35010b93 GIT binary patch literal 20456 zcmbT-c{~;0{{Vj5wXfGs_I=-X*StdZ<=P|b4I!bl*egUK6iJpAOC+I0)<}{HAw-Fg zrIfuW;&=KzzVFZP|DS)Jk309y%$b??Z zbSv&x7Ip$!k=KbLj?VC!Bq41W>jEC0I3BJ1gJ9Lj9gzRP6@&+0^H4u_gE>Xz8s@4) zaKu?w&Mw~_L>V6x9gp^aU%jO>LN=a|`Sk+_Kf4!{Uvm!9cv+T3akp~A zWN#RY8^bU6dILVvWTl_U2X1(}I`WwMfZa4c0~h53U*k#M3wM2>{uwJtu+IlxHx1UY zuK0kO*U7g#*o)S-x-MKlE5@3<*+x>!#5@>6Uk|9^{yVd-+I$@7|9x7;PeAY{pqIy z*$vblPRZBINE@-mHD09d-!V&3WH;pzH)-U>2Q8`cNX9pwI!ef;;>k*7z$aF6 zVdPxdDLnz?T61VPAJUw=xr!SpuduqxfppSt1Uxc5ccqX8IsCAPfeC5jx0ym2fAu;#m2r24OO8()wC7p@Ec^E}12^9-`Y zUFG^Ta)aT;=P9J#g#&8;RJ>=B{uw#P?6W$F^eRzt|L6aOF_~xC+sOI2a_rXP9QG_? z24iX%rejqhg946MYa3Tt;bOg3PW)#+c>mz1LxH#$?6}2i>nX{AZVW@dNthCxp{pOu z{H6+iKOs_0MiZ`*wkRCJbl~wX%2|h2ec)L$oD!k|P_m?Y_pD|=$i)jA&bk>vboc4l zb*cj}8u_)IRAmCUjk8)FI-3E%@cZk&f6c)<_2LRektHxMhGlXYTf?v9OEt1HHb5q+ zHR6X6fiYq_(TKqg3~H`?{5fO?Uk@s%9VoB|g&QG1|9CqsE3!OlTm4Bq!$r(~&%BqghyTC91o6in4yTFP?JGooHK@j5$)-qsog*Vp*1U}cf z!tCzpL|$h%a2Fjp*|q8hYP+uwy|{P?D6BIyH~HLw{>5GgR_+e5*8PtzD0{%-r0&$C zS`VliNJx3C;t7Yn!V{w_JYin@?&%!?FKBYGU{g->0zdM5!2Z(<1fMrOaFF+gK<+Q% zk%zot+p)lI59Oh&w^Tp>c<2qcWs)at#=OC)(!N=7*BdB<)fEbaeW1fkM4}vg!2VWX zxu2&G{Jz0e_BGN6rl~r*EwX%IH|?RyoeCcyvkcl#-R=W99;d9b-}yklm866JtPiwZ zv3hd-rw>#Z{#-MnB7s?hntc)*3C{Muco->w`d3@Nk+w7mUgXdp`lv#JMf*L7M)^!p z?JKL2A^M(;?Ixo+36A*&`AFN6z%_y+^qM0H-1Qm5{ai_K>-~?~at{*x?4aWh@j>}+ zA{w6iqIg?+7xnlE30_ZsU9k#6c`V=mK~69UxKpKT^Fm0l$+pVogYx1!Ec{tHj`}NQ z&JZ!21j`BJHc2N*(8JgkmUD^(v_HimA{ngXAY>^bwYlX*ov9N=JFx3~(66MTvdh%v*QVRwAreFa|^C>)oQ%!oY*7V`3M zIZ|E0n^kk#FU<|yd?#)yoIeDef7vC40^Py;#qaSgBM(@rh)zys^91IZ$MW|EJV9A` zt<5&k3$)W$o1bxb!;fNP*GJyoU^9|*i|ej8JgC+Q=b7~eyUP`qmIQnte3{1Z00~bDJdz#2xK^=edyJvj=H974=uyus!#&KoSJ6{)p^A{qdaWmy$Tt-yDw! z&>NtB@i=2bSoSQco6oOU`p1(X#5m|Y>jhMIeTO#pCzGJA`gJGECDcDzQ%o^7Dc=Ig*lcv$ z68X^5lEVsFLr7)&=X(~(s(;HvmnwGt{mxM0KI1=)ud6TqTi;7qx7ZX}D>qO0w_P)i zxejCGBSL|;5mHIgqh$ZT<)^dh0eQ6STgbom%*K$0^pPAoCkT4~T;{3NK?)vN712h{ z&%td?Vlq>!C{k)%mGYkrhwJS`ki;UBa$#hs%~k4u=H9Ug7DDcw6d4mlPF{Xv_0Rv} zSSf<^dynu0J7<;A9VfPqG0rk%EA2AD>VjKFFogib@a6T0~J13j((`DZPLAhAXtwQgh#yVu?I zXsZuE;wi;6vA-s8hh;c}T+s|HGqyu;hUTEFtel~uYXOO0>{=OkETQ}&7vHH*mawmD zm2@`E3RZckADc>9gUR8o=?~S`u%QslD=uvVtYUdoU#{3d6!)fV)shX!eqG%64YP$s z3(XfJNwyG$dsCZ^@=MOI%m<8rY@v$v0|O5S5tt&*o&PLFgp;q&e`wSuLcOb9|4$Pl z?5eHJPS_J6OY4DSqdO7BCri_FUn1OJcz#9eC=uvZ@{-NNh!8fGP$hGc2x}3V<&&Jxzp3o11mHr-*RD5x*dZwoh81eg5zS`kr$im-2BUBnY0#x*v+RA6x!c zDuf8QK6~5pU?TWW>3vc}?dv$N5#Aa|gja;i!A=20X#RPrd()o?nzWI6rG7-9<(t#@ zM8~6JeTZG`FcFORm?zeKh)^IuUC{4Mgo?Bk%|?1vNZQ$8??F2vyk8*61|#QIZT`?Bzt(kLts)YEp@*!|87byGqE3nAcgXLr zwI#xz5_Z=BQ?cELsy`DOqN( z9{G#MxIGa$7-&I;?ubEOz4De6@=lGueb{nAiWa7RGeZ6s$9$Vr&+XJm*!;4{pu&@CtiGY^*sJj&%!@y^7bo`vqm_=e zTFd>53DYIFCePAB)3K^iF4`dJY}Ad6bw^WrsvK52UM9cp6`Pd^^V^vxL3 z`}1;mIu3v_$t-y!-2|%QuDwzvnu5ZQ(gG7|GjO>b_R{X585o4mZa%a#hr{NY^JDYo za4z$#w)sg5=;wP=ePPoA;sYL$bDgjRyImQjxG$E_?s`~*U~2`6heA8wJhB3c4@au= zS*+mzcij-Tzcu(@G_x+Lw+4B=D{fSN+NG;WG za^10k>BiePl$CA4{gX+#hPN%qW(w^Tr`f_v@k7SiI$JnB8nd1A(H07l?~b0^wS~e- z$rH9L=ssdh$e&q|2;VeH)~?8*d}2*k^FxaWNs=l0>PAF3_&DC~4C+@2H|B?5qIfZu zSpT7c;<3k$k#+#Z@spcO^JT|SoH=q zs|@$A5#hHOmv?;*5wcYWG%sIAdGK3!hDI(Cbf0I=KSo}ZngG@tX#WNoQdUT}r_7DP z$ot7ZD^4MmdIar*kdtk#_pFh?Ln?&vGqCKY-hx|JnQ7wRs`*nTa zQ*;~;a&^YMknw3M%9O~@{B54)SwwK{vCgqY_C@~vyM>OsXPw<5sbl|^J(y}y9dCkVBk~Ox(Nue zr<0Y5n!)#+@=|IIb2w>j7MCw)0l9(R!s|4az`HM1^yeE(D0uU;y&=*HKA->d$BW$> zUTnIn?c`ZQ_r;8r7`zRdr=s90i?)IO_m33uqc))0)fF!zV+;AQQ=k-P3ky7#bD_!> zss&{^+`iaCoCw|g13Ds%h1Qi)qqrEYlRloQLqlz$CIyhBl*(^d}du}1Mw zRRj6QQ9Q^=PvPoqzu>*FRqm46!PZRL$>&i2QjzaxK2L=9>RFO0(&R~w8ive0NSj1~{AkwKlaGGK zu0k=zjud>tpm_qVFP^|b*NV3HrDHSocPxr4)7uR4$P<5X_=9MFJwYo8=TO~OWEpNM zL-l?=MLqv5YR|syma^XwsNXl_vq>Sn=?``tJcH^@t#)%i@(uvif+C>>V*9JiP`?iaUyi?Hl_$6Gmq*0xg16WbDPpQ z3M5~(!fH$?T3&Tj>Td`Uf*na4@yNj*7a3+G&EAEao9O&WEtN8;BO9Ol9PA55=a=xv z>KJlxO{$RQLMb?qXtpZ(We^(2jC7a9A)AVHb##z1GXk8O zXdJ?Ou!OS#saoXE6@xsUiXj)Sn+LKZZJiQMQX&6~r>WRV=L0p2w}JzAIs;L4Bps4c4lh2NFbjdXPoXu(f%il_t2Y5fN`%G5wp z|HwA0xEj2=Yai5bT@{Q?NR{c7s*tTme)&j%3h1`-2kmqz!_FUjhd)Bf@YnJ!+r9-w z=dH(Bpl!yASXn&;N z5r?tlq2>5oG3bn`*={lrg*=}>N@rGuK~J>cMv1==X!@D%0h0iHR%yL`?*kwFK4{wa zet-u&uNypS+28`L@calZV@`z(1VWr71XBf$>O130K3cZys$!hzY!xNEf` z4pKRX;?nlxz%cv>a}wYnJCsP1i{T)eM8nf&fP)W#ig8;69M~r$eBY0(*Ew?PhCU9o z_AW&J)x$yKlUFx_kpeMG_h)r+@O3)q;Sr>i)BM|29UPn$p|MOr?iW3LmLJI(ba0|k z8y(L=<2Mha$~PWO8l+lgjAxw|4t_3$GX@~Ji_=#Hkz&Voy2dndpc9>8m5Y?|8#8r7 z+FKt@5JE0UYWB@*;K1k7_Mb;c{-3jx=aAoN1f1=Vth48I#F1BMi2=XWao}r(d)|+< zo8^&T?z^dYB;bZ(3|!m{bIr! zYmxNSt0(f2TCa9ViOB!r*z^qFFOV+7gbsYZd0Bi4v*ElQK4|<0d$5oa_^E&vro>*a zglpk}vgppSf*NkfJ9Ae;E`A?;jn=rQ=^zZ}bPLnFIYmK>I?%uHsu+AcKKjY-ojBOk zw1_!>mVlnn-LkHBNwARH*}WVk1uDYVb93?1FwJFxr@kx=8wMkX?SDu^Sl_u5r!{4O z&5BQp{HP45k~4jly(|OvS!3EWH)X&>WISa2jtofWew56*E(7Zi4t6-6mw_|y)vg$O z%RqN=!+CEF8EEv8ExkZ41DbM;#V21#gSZVJuL!t{OiC_f2t)XB9o#Zb7{-}q+c$p+!aIqW-03L+@Kdh2c(ivPVEja)e-l4w z&~fkGzsCnJt$SE`3wR-ydyy;kDi4^Le`1=r%ng^%N*(!{&IKC@4j*E2IN_U^jJQWB z2N)I)yoDBaNU*V&>HUcAgAON0XHnsySpHNaGs-vbqADuxgtLNf(Z*4^Q!FqU_&Y77 z8q`D7pv$zgeu$eIjSrtST{wz^4y;rx@hufxT;RD)G^7HPqB)-F4ocX%$*`kLq=dBY z%T2mpDZqa5WEs_E3Q#@EXLwVW0tTqLeu=J%xSI^stE|WZ-)~Gg?rO49Y#KE(gk!!8ylMGf_fhU?$i4h+dHGA08D5pn#6?@n>>8_VwOZPC3&7iT|JA0T_tHJE6!aYo1;j-wl zygiJAyu112^*wCrM^ENZ_8xZbi7#3F)je$JOKC%R<{s8V=6Lx|`X2V-k&gq_K^uA9Dkzsi?(#UB-rDYqT2Fwz`SGT2qC9PrgbzKB{2x?tq4Jw+fiN zj1elARsrVwG3U-FD1*Ai;Hq$~60qdo;<9m3gy!)P4&9&fa4EF$`u%h{P?DJPPnM8{ z_V+i{-Ze{uNpk@Hv4a#efL~F-vIKA?T;s#8h(p|?j@o2BG58~L-#T$c1kDfI9dmvp z40ST~93OB(aMtvZ;Xw!iA?_HhklPpxZ)|e~eH)PGfe2 zPIk~x8jiXhgol^KCP_saY*4DhsOKui3XOjX&XOHw2FXJ18I@*6IIU`_L~p?Wnt4ta z_37#1cxCaatG{TWXu!{ao}L!&M4T7=fzdz{za`Z~7B%d{BG#tJsp0#_HNnG4I1rSP z$0-Zqpn!fY`b-xU;2$)$cEnHt^$m8lB?BthCil1|!axO0HPzoB&2qZ=C~ z_{PS{H8N8|OJvOT7)DB1kVsg|qNfDf=h}BxX(@r{d8S445vOE|+XH%kfF8S}9p`Zl5q>xS%jKCZMQf#7e<>iig-Le=nQ45Qs~%Zc_7azi6rS=Oi$?aB2M72d z`Bp;H&5(JQT;)`daJ^QPADKMa>qd0*sp!kMkax)3BCa5h1T70DA}vgLPoF_D?z{^>iu~+xmDwBlUmR=LEWebRYR9l2 zxb)_WYgoKes{!|A8YpRx8R*q!hqSOmMb{H}q472II;XWD5Ni$}QT!kbsdCLiZaX5- ztDvSni|P-1&_FDaLJaO1_pVg2i^1Og>EEr(qJS4#zCl|o3RXGNr#iJnfj-tG>u|jY zL@RxX{w^T`x~YR33TeVH(*5*V`$r*2X;e3hFBF86=7)Sd4hw+v3zyb(m3<)kr;m$+ zfgc(Q_t$-2mFyw-YAghUr~mK(81VoXfF0^WK>giYSj_ zJGmS%LpBxk{VO{JZWNe0n6tx8!_Kh}^>`4gmZwut#RKc5>6*1BHlSLJ*Bcpy1H1B#u(xQkJ_AhQrlY6j2p8pMH9E#HbTi3$ac@}6ly&cR{$O0uS!l^cXEKvD{ z#{Q`R3-nL0+K(?fzRU}r|fVJyl5kH)h;#4EFaNQBKf1~V2&iPM=1KEeWN zp@x5KF0+7Sza0K_0}E{a6uaDswy*lVq+5WA6%4p+onxg~;mZkg*IS0Hu#=}nsqD!L zHMN2J7b01~>843`+BH_ttL6W=Si=gxpY4iFykUjNIPz$tIaWC7+VYixj18s=_k8Ji z*+62jX~j&14JcZe2cDR*!E?%_VOuXY7yzecrc-S2_K!yC$|W|~rg(MtcOe^W_A&X$ zJ!FH+{aDqd=WH;e`rr`nFdGaS|&l5jt2iia#?QU4)(Jc#2n-&cF!;m?2^lFZBd%_4Fl=K3)KYzf3<|C+~n#6<2`ND{eX*@7DHL^UK!vk^b=Q{Qc4_b)s|GX$+){EY|KhUhfO9XuFm&%G?#MdzU# zZ1>36;qKtCdIBXo2=bYo^}?~k{puJ_e;Ri5e#_Y*dpdSdmyx5EVqk}5s>YcaMs}!n zT**mcW`|d)CA#}q+5Z>Ey&>(xKbntWY0F(kGYk*0ubggE)1*P{V0IeSEwMT5Zw<}* zrIvN9UB5u%PU2tedMfKZzjjKf=VzIxNuh$ zkt9PII~YBBDD!Nd0|J*)KWD6>ai9F>k&oS65K_ZcKy`u}bcf%#=&y1Ez1O&Vlm!o* zv-@rweVzwGa4m&hH+Udk;N!Kr93BYEDYJ-(=7Fb7&=tX z4TfZfhpb(=Ay}!}q4ozCB$ck*ZHwlD)ptj)hf;FE?R%0EwyMv0cKEryd7hn|9nFV6(`b&z!?*4IG}F}R`YAJc{k+Tu3HD!7 zEO^);C2IAO=v`K*|4b9~K#LW|wrRup@>szBg^ldP9#n_x66Bc!1nrFGTWxa;Ku!53QDi>@Bpy6H?~Lx3nJmbr%$VpwrByd&*ItPZN}|A4fUC zFhU8Ulz}RO)fw5QpH#s96y$elaB$Sg&0v-n2hVmN>G{jy0C$bEDH7d}2lfh$^P>BItM&>B zLk}FdT5Q=b2cvlbH#0S%7#wICo=xgc!9h=7?I(&GIG`L)j;p_mgXHDw=HF}4?-dk_ zOWTnRv7TRg(0VbN_bWzm&>kU!+n&L}uZyh`qf0p8W76I<{*CT`FXvn`rlf}5odMbp zjMQLw&;2ni7d6yQynX*pkQ%0HN;(assDW2zsK-i)8se!i19C0u|HW~0?veU`^*E)r zdp!F#0R{x-T_^nT(N%~m1; zMj@fHcs3E(dST2va8ej1f8JSVDHj4SW!49?zJhRT?m*3$gaCxTmrf8`<_EIUY}ShB zd_Y`LdK_HJ3lTEH13A}tVEJ9H*vo8mpIFCStx?Pc`jyQWhuS&8?@8>1fFB%S^>vXa zUz!6jhv$vdC)lC=$IuUhS9qxVCP6){fQRzhc6+y5Y=GOGz2q%`>Q1@gPDB|iR65^e zS0k{(xpEmPjTRP|)+}8(%FY7EOs#e4&oYCe^EnT4G_Rkub9zB7g9$8@BezbgFafpu z`Rwc&Mkw?Rps+1rgf$nncHa|>&~!oeZMGF7&_?l56sw?dc6JAyyD%f%xH(FHgPRe2 zJ9P34xfr3J)0Dx1pApiH+n7cr7}0y4e0SE>8DWtAY{nlbeim;mZ&oP40*$+d7R7^n6bfoyKC=-w|ol!Rc zH1GC}wer3P6J*7++$xG<0w=HDEvoBG@cP!nKTVIAz|@G+`B5(uc%Ode`DvaB{Lc(H zsgg6p)pGUM6}-&A%0b3#q{0lDIH{+ zSvalLgnV>jE+`u7?i8Cy)*Kb(Kb&ds2T|D>tYcdN69zJI7 zp1}gy#>Z<;q zHm=_uMccg;&&5U?TENW#w`MVNxSg?~GPs;sG ztnBv{*Y;-=z+7>bZ;gcx+AXiF$ci(A|Jjop2fwp|oO`Skj~E_e1+FG5zQx1-RqHuZ zRR3evZk|rqz(cp>e8Iz;c<|bIGFvZzhd|1h-%;n-Kw!-%Ms<)C3^{0xzw@y|Vg2@Q zLOU~z;+Coztx$i`&r>M>!~knD{_(PNIbXuct1 z1J+*u#g@24gn+bzsR!mB59ip%CIqW+FFyaqu4;F`4C&g!q)R?;mvwGpiOR|)w?F*C z7P`w8V`zS1E6s1r?JRy`eq0|Hg^D(?&coO1-ZO1r-wk-y9n;palKbX25AR#Y94Aay zJ6hH-o{6zJncy`{d&f=JN@fka)OqLpFz29991Jjmsc^R zwud#;v8&ko##_C!CsxrjK6}CVkX4M~Lg*%0=qg4f`;@jYd==~eGyT3adKLTi_wvWM zq*aX6l0>|jwTiKLw@jrLtzrpn8k&=jR6kUF|3E|S_p~+a;{Yy0wr~y8)2I3Psd5c_&`(^UM92M_p6qV^>oqKGN{Xp; zbPcZVhV+5Zz;2U&Hopw111-TfsU|@KCD1w9ix5p zHh@AN4Pif$(pELrvD13(e~w}6nAwc#M6}sD7U;0)BV@OZ75>R2yVa(W%BICYdsFm@g5zf!tlmVm~YxA~o=FRf#vO5bm_ zrLSXi~*X$s1)>sX9e zawJ9dI<`khv+8`bj*YY9L|y9Fv5?!3G+P_d@eHj#*h1r33W}>&8QRyesJZ(0n;q-e zltQ~}+lzIqW?O3VaQ8ZP`?)vo{Ofg$zNMkmq8Dw4SByFH?K<{4P4IQYz&a**+Yjq{ zw~mRZ^mjB4uVbmDTe;cq*DNchZt-ewysHB1Wth~khUeW?S+i2#^4LUe*@vptO96hw2 zYWvF@S>smDYqy+EWDeTifX0<$*7j&CxH&` zTe8?V)zL!OX|1R>0b1Z#$}?$>p#f63(xX8%k81fr{FJ;jHL&nz3(Lz zEKY5d^P=(F9Q%xC!x2ibm#;Zn+e86=O&NFRSt;Q6S*Bf@qvR0vA@ty#S7gBB8@nc; zKnB+$9gMQ%_AqDZUGqwsU5uo1ss8EG4wg5&>pSsz8`HKg4xyO%gIPMphB3_j#_)V5 zTish*nE2mBfnClm?A5G;yB}c_Gv=0iUK{cY^UFP+kbmzd=51NO+xue!%iW*h_g8lV zbHm5ae#%_O659@6kX&BFs`L#Tt~jk>>0gs8K0RN>@|{ZEdURH?Z!EXt#h$ESyYY{m zEA>_|c@l+%fBG_(UU|KdsdfpwM@*q!{qh6L-n*n%C-wt#-Rx5)5B-kCJ@)mC=~~22 zgl#CZ2_YhsaAdJ zvw$(W9_}YGEnxMp>y2xl&0_&k3xfT4hz|i zn6hr1!&cfZCT(Az!zkDuQ}2tJ!xrPbr49V%uxn~HJRA0NSd+9O<&g0lrfl-zN`&ql zrbf(iJg+u~+0=8gW+}~KDsRL$TomW9$xW)3v*>dmbw*>G>Kta+=H54`HHW#l`8Y-z z%wbjc<$gz)%wb)=OP335=P;J=*1lTTIn3y7y=a#29M?U&)G2ODY!Z(z>W zv5RvUZ{QI9t3+t=Qht_oPT6;g1hIizF`#=iQzeH zD7yD#{nQ*bUU$@kYiSN6`yQR5|7Q+MW3FKprJl!Fm3WEI*ygdjwi~^i0`r(8-3PXP z((~A}k2FaGs`Hpr!kC*n0WA+>J7#1$kFDK!s<=y>#~xOkx$5gSj|l`!ZJj`2aTTAO$E?$cK_4u#G^z*x@9=qJ3+q{A5@$l}k0S{D< zn|rk9j-q-TXRGYOf$A~7tM{pwC^OhSI1620R(ANZ*t_kxj(F`6gsq%FJQz=-DmdX;Xs@{_||QF$K`T{6%$i?b5x+r?6f zMl|J(|6)eNLVG8Mw=u!htpil&|6pBL3;0$J{Kg)tHK?13ZDA3aKl%lEHZiNXJ1-K& zeqqIjl0A}*e_|$5yMBsM8`#s!0k5N;uVaeGpD+0dtz+3K{Zmt@KD|8p$NK~Q8b%!! zdQ~%b6&q{@%a3ZS7`{a1;e+uNteWBrWzLlqZ1KK;SG?T{#?R$q6Ue`U>1|??x#P=N zr${yD@vF<2s`*>uHRoll(3FMjkkm4kFOoZ`MY)Xe7uuhFwy=beReUOZFtUVYZj+Y# z-z;HNGqxInFPE_0ovV*ZkQu#MJYUAzN4z)-i{sI@LGs8UCoKA5UsCyrT8#Dm*L_o}pRworMLGO>G=On1e$!S6!{*8hoqZYg{C=?8NsRdRp zMhEov_6!{Uj>HpjdAaP~!f2IDh*cr$qR5#yLX$WXp;Pg>H4Z(Fy0{Wx@f}~9u_poUuYXiFArauBZ%El0KLTinC!Jb5LI8PpPyW6@0z6Uv z&MJJA0G8>@J=@4jojrXgf(bAl>D7{ryr-5F9*P_wR?4#=9mqoS;?VCdez9P98AJf~ zldh>NXghD@h&T#lfe$B%Jb(c1hhy%p`xDUf+X?1RXulNeKGM&83BY}PcE0E^0WA7% zJxxUI;(QuW8Q?>JL(PUnqBjAuC~1D_dlA6YJgHa-wZDU$PE5vw06h&ojI!QonJrBLRw@eSG;7ov+qYPls%gSrK!}5A6wX8W)1{IRU21cBU@c5rE)(u5p)0 zfF`PrxhSN<-wD~jwgk|*SZR6z`Rxmij2+2+nt|u84FUash@F!GvYwmb>_=;~UGv`_ zu}BNikrpxJ>MQDiJ}Uym+w?2MB7gMQOei5=_=a1~SrQ=k)SlUGWW>qXN)KdDwzjn( zviPA2)wBgVj+-2n<;dkFS&~8IMxJVdoTpi7+yO0U z;Co$X*`@_AS?Pxj6}6$HUucgaTN?zJLS^2rYJ)1-z~(J89SD-`_R_qe1LjFdx8%O) zfbU&dyT2%3jQa)-a30l#*(2P~2g-Eemm_~h!e?CwU!%%B#i<9+9G?~~n&^Sb&Gas@ zQ+gnw63FdxTMv3JW=S@6>%r*u@u1ixJqQaW<@B=Y!#KrWXpugwn$TGsZ_)t$^}+Snl%468KHP7M zR#%}XfYu#N?P)#&OxQ#XAC)CQNrKYB32hYTEm!xz`u)370ci}y-S{uh z*)IVEV0jWA^5g^o+9W5fI-&`n8C8FYHIV>RCsp#`97p^;3QM^RQ9^1?%fTH!5 zVdG5#{Q2A6yrnNHN&;*y9i%I(L3!d-W-#?*6!+#T{8Uf=-~O9e z)<^)Xzv0zpGXW}TC>yq03GklN>vdxr%6~aLlG@Ku|1GA;i|HT$`SfmK+H;hT76$B( zbrPWQ?qqt%Pf#D&8Mx`Q2QB^j$X+_8f{VZ7$MKN zE}AW&qb-Hkx?n9%RFP=@Mp&4#ho zCsdEgFI6pEK=t@;pR=q5s>c_@U8;|vdOR$Bp|)|KAaryKv&9UHq4(~Oke^PLhDJS2 z&tvGl46)ISr2S`=LEbOtm9?oVM7-3Z+2~UPUpsx@b#4tfcsK7zJc_$JmLcEc=(HeG zAUOFPL0-u zqWKBogl=8PHHf^P#Ek0h{?JJW8$I|XVjhunQ4e_N{Qur*Ma$Kh40%vIF+bx~vF6hU z&Pm#U8;1JuVuO$03e{tE!KzD@s2)?^%v0P(_1Nap)%Ft|`rz9x!EF9fA5x0ynyXh) zJaLE54&eys{~4XiX1ECuBN&%TC5`f8X!mCiO#*ztxVPDj2r!iX-pd^2eU`UhrUX$P z`u$hDrq-JPyJfL%5-5(=QX~5_LI}{8GVIQDk^uAqF}Ld@2#^%ImtGV@K+mi9j?JAV zfXc4*=(TtP%q3{P=tv~MPuiiZ*d!FU^}>9^$pkp~Q~SxtMHF{U0s-ll2(U0?KS7;B z09x~s>*h#m{Z6qS^?!6yjCx^ws`NaHgO03BK_r`D#X$Z!w45c&RRHO*L8)Ji zj$1zD)_Db_uE>&1OB?~V>c|vqk!*eMV&-De_-6Dxbv*LzvRSPNQv9xO;0tu#nt!}u zJchK~_v$PU(uFpv_$A6yLw>pkqmbeE_hOZi(*i{sU!zdH`^8gm7db?%cJxPfWk~Hw z|69H_9kUdP&g;G9>UyNdtYGwcWKkbShYixq)hSRAc__=oc?Dg!`96)%PNe$}ulzhD z_GY2u1oFkDk9}51(f9KAq>*CW__xPD=u}(8s4Q7^I47|aFB!i`ulq`C*;s`M~|;`sd5nSyfqq;CkLujYP~bJi9DE+;(ogl5gF9l9+Hsn4Ksq&Z;m=Q5;Ma@br7w{Mu6GB@eqyq4 zGX#6$^hyTO>Z;JwF&TJL*=yBVA`MidV<&jFq@am>f+|W+3cmMkFtmh8Lg*K#Cc!)j zI5>DUK)*>GW|JAISKf(17|kTP);CcYbJa9tBNK)G;rIm~9ubIC?m6*JLm2urs>0&y z1VJPte#-_g02Laqf2m*KhngNiuZygF&>rMb^0tl#@S~fr-ksuxmogkZ^+Yb{R(n>; zWWtH+N9mnyM-K378{)8zVTY()KIzORG|y1THz3Q52OCwnY?iG1a`q{hk$S<@%Zb_KCOOyZ@mQ1l#IiaBfF6GLX$I~Z)F@P-LaKYuD!aF_`; zF8x16+-W#fYa0jfGM?5lW@hhg@0}(^=P2h?#6Q_WRLZb+hUBEkP)Tx-QufK1p@W=M z8Wf6@29;>gF}~C;La9WAA~KYN-lq@m^?v$Y*P5>NW!=xUuKRyK<$W@gJWh?CZ@5q@ z2bGE?W!-iIXDcM7EcHD_w4eJ4%kf~Q@%UQQ!>PHL?no5uOF8B;bQ z8cC1u+~Nn)SUFjF)oK%sgOc9_<~B4u=bTWwu$;!vHQ&X2LmD<0+3CKsXsAlBms+Mx zV|caJzsyV;pMY*^j z_3sbT0g&{w~Ib&TE;%p_8J&0~&ZFgJWQ(#qR-SQz8JgAzKCURs(y^A_2 zxa-+#Y$L9c(b+$n*k}5O>8B0~l6Hv;ONpr~{#=_vjMyct@+0y=U7Cq&qe@Rz5ha#~9=}8MebFOdNVL~At3OLL(LG;u zoLCq?k(@;I+gq~!2=RaGaqp;euIj=Fb8E+2ntR#h=9L3dDFdHhn;X;|*86A8pt<=W zbG0=eCd}(*#d($Pmckm-?ybzu9i%BRMKyt28|QY@vpL zh^T1y3+fQN0PUPY2_u45DOBdH_hkcqfmBUfEme6 z2fFq*tm9*`w0VHgo#d&#ER)}EPp>P9~+^RY*##W#-RMfs}Ps?Te97&Hj^U4`Vc zOB}@vL;g=r$JE8$GvQ%Z9#iTg!9yeO8R9Fr=)K;-_BzDHu1y`C#-?0+&GU}9-l+?5 zV4u|{TU`VwSEen<=V15`wPYD-4(70Yn;Ko&sDI#~gH! z7Twk>KBJ9S+99?LLt5};XWsS~XrUtEdD+3gH4$tSG}$nu0i)1P3AuF|sNIv2|MZGF z3NN;^bJf(*>hr{VT9g{thvQ?7hg6Y#vWlDNtBS;LbDfuVs~}%w`;wLTHd$dKBr?hDGZvtPv{WKfmklknB@H5`StzIVZ`*));hE8pK`Q2zKu5JJF z=(dyznaB4xeOE{Rd(z61yRAqW67JCmKED#fs&;Lm-z|{dIk}(BSmOMGj9B+V0rKud zyjPZ3jrRVz(qkK}!I6*CxteW*k0F&4^}V*Z7JA*Z?q@rkt~UERGsqszQcKTI7OaKa zz$Lwwj&-oNeN?a~mh1yyI5ztiwgVdP?HVY*?0`oH8UEehHUQBjf%orizyZhV@*uN~ z=v%2IKKDc>OL1dX^FlfP-M^=u4e1Mh(XHMb9Ksg4D@oBb;o!!=+V*|)Qw;u_jiIj zH<>|Vw9k(EGYmG<9aBD6$@#Q%Mf1xU(y!C+vaDrDp82stZzBWa1JXUK-Y|ImHC7zj z&A=ylAL~*-18a`imT$uhM0SFf$Z-byt&U#a_?>}2OAm>XLWE`>i!_rJLf>qqch=!@6(Lp#a$`2D3DI_^D!p1=2=hw~d!;mm zIN2+ixuE@})zoFp&=ulsj!W=t zV*W(@h84s%ea69oXp$`=4P_x*nnZVY61AwTY2idY_Xd76alikuz##6!P0 z9Z4ewZz+6moLJozDsh5X5Iedvi+Go@TyT=;@l@8~hnI&G$A26j*;Z=u!-juN-DZVOK!`d_H^ sh$9LdPppa|I`;N2{F8X&n=$*xd*7=~Qa?oGPHg)3t)DvApigX literal 0 HcmV?d00001 diff --git a/tests/test_text2model/test_error_message.py b/tests/test_text2model/test_error_message.py new file mode 100644 index 00000000..f21705bf --- /dev/null +++ b/tests/test_text2model/test_error_message.py @@ -0,0 +1,57 @@ +import os +import shutil + +import pytest + +from biomass import Text2Model +from biomass.construction.reaction_rules import DetectionError + + +def test_preprocessing(): + for i in ["1", "2"]: + if os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"typo_{i}", + ) + ): + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"typo_{i}", + ) + ) + + +def test_typo_detection(): + for i in ["1", "2"]: + typo = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"typo_{i}.txt", + ), + ) + with pytest.raises(DetectionError) as e: + typo.convert() + assert "Maybe: 'binds'." in str(e.value) + + +def test_cleanup(): + for i in ["1", "2"]: + assert os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"typo_{i}", + ) + ) + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"typo_{i}", + ) + ) diff --git a/tests/test_text2model/test_fos_model_from_txt.py b/tests/test_text2model/test_fos_model_from_txt.py new file mode 100644 index 00000000..3cb5de40 --- /dev/null +++ b/tests/test_text2model/test_fos_model_from_txt.py @@ -0,0 +1,107 @@ +import os +import shutil + +import numpy as np + +from biomass import Model, Text2Model, optimize, run_simulation + +from .C import BOUNDS, CONDITIONS, EXPERIMENTAL_DATA, INDENT, NORMALIZATION, PPMEK + +PATH_TO_MODEL = os.path.join( + os.path.dirname(__file__), + "text_files", + "fos_model", +) + + +def test_preprocessing(): + if os.path.isdir(PATH_TO_MODEL): + shutil.rmtree(PATH_TO_MODEL) + + +def test_model_construction(): + cfos = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + "fos_model.txt", + ) + ) + cfos.convert() + assert os.path.isdir(PATH_TO_MODEL) + + with open(os.path.join(PATH_TO_MODEL, "ode.py"), mode="r") as f1: + lines = f1.readlines() + for line_num, line in enumerate(lines): + if line.startswith(f"{INDENT}def diffeq(self, t, y, *x):"): + lines[line_num - 1] = PPMEK + elif line.startswith(f"{2*INDENT}dydt = [0] * V.NUM"): + lines[line_num] += CONDITIONS + break + else: + assert False + with open(os.path.join(PATH_TO_MODEL, "ode.py"), mode="w") as f1: + f1.writelines(lines) + + with open(os.path.join(PATH_TO_MODEL, "search_param.py"), mode="r") as f2: + lines = f2.readlines() + for line_num, line in enumerate(lines): + if line.startswith(f"{2*INDENT}search_rgn = convert_scale("): + lines[line_num - 1] = BOUNDS + break + else: + assert False + with open(os.path.join(PATH_TO_MODEL, "search_param.py"), mode="w") as f2: + f2.writelines(lines) + + with open(os.path.join(PATH_TO_MODEL, "observable.py"), mode="r") as f3: + lines = f3.readlines() + for line_num, line in enumerate(lines): + if line.startswith(f"{2*INDENT}self.normalization: dict = {{}}"): + lines[line_num] += NORMALIZATION + elif "y0 = get_steady_state(self.diffeq, y0, tuple(x))" in line: + lines[line_num] = line.replace( + "y0 = get_steady_state(self.diffeq, y0, tuple(x))", + "y0 = get_steady_state(self.diffeq, y0, tuple(x), integrator='vode')", + ) + elif "sol = solve_ode(self.diffeq, y0, self.t, tuple(x))" in line: + lines[line_num] = line.replace( + "sol = solve_ode(self.diffeq, y0, self.t, tuple(x))", + "sol = solve_ode(self.diffeq, y0, self.t, tuple(x), method='BDF')", + ) + elif line.startswith(f"{INDENT}def set_data(self)"): + lines = lines[: line_num + 1] + lines[line_num] = EXPERIMENTAL_DATA + break + else: + assert False + with open(os.path.join(PATH_TO_MODEL, "observable.py"), mode="w") as f3: + f3.writelines(lines) + + +def test_parameter_estimation(workers: int = -1): + model = Model("tests.test_text2model.text_files.fos_model").create() + assert len(model.species) == 36 + assert len(model.observables) == 8 + assert len(model.problem.bounds) == 75 + + param_idx = 1 + optimize(model, param_idx, optimizer_options={"maxiter": 10, "workers": workers}) + assert os.path.isdir(os.path.join(model.path, "out", f"{param_idx}")) + assert run_simulation(model, viz_type=str(param_idx)) is None + simulation_results = np.load(os.path.join(model.path, "simulation_data", "simulations_1.npy")) + assert np.isfinite(simulation_results).all() + assert len( + [ + file + for file in os.listdir(os.path.join(model.path, "figure", "simulation", "1")) + if file.endswith(".pdf") or file.endswith(".png") + ] + ) == len(model.observables), ", ".join( + os.listdir(os.path.join(model.path, "figure", "simulation", "1")) + ) + + +def test_cleanup(): + assert os.path.isdir(PATH_TO_MODEL) + shutil.rmtree(PATH_TO_MODEL) diff --git a/tests/test_text2model/test_model_construction.py b/tests/test_text2model/test_model_construction.py new file mode 100644 index 00000000..ab79ded1 --- /dev/null +++ b/tests/test_text2model/test_model_construction.py @@ -0,0 +1,184 @@ +import os +import shutil +from typing import Callable + +import numpy as np +import pytest + +from biomass import Model, Text2Model, run_simulation +from biomass.construction.thermodynamic_restrictions import DuplicateError + + +def test_preprocessing(): + for model in [ + "michaelis_menten", + "Kholodenko1999", + "michaelis_menten_jl", + "Kholodenko1999_jl", + ]: + if os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + model, + ) + ): + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + model, + ) + ) + + +def test_text2model(): + for model in ["michaelis_menten", "Kholodenko1999"]: + for lang in ["python", "julia"]: + if model == "michaelis_menten": + mm_kinetics = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"{model}.txt", + ), + similarity_threshold=0.7, + lang=lang, + ) + mm_kinetics.register_word({"dissociate": ["releases"]}) + mm_kinetics.convert() + elif model == "Kholodenko1999": + mapk_cascade = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"{model}.txt", + ), + lang=lang, + ) + mapk_cascade.convert() + # test thermodynamic restrictions + desired = [ + {"10", "11", "12", "9"}, + {"15", "17", "18", "21"}, + {"18", "19", "20", "22"}, + {"18", "19", "20", "22"}, + {"12", "17", "19", "24"}, + {"15", "20", "23", "24"}, + {"12", "21", "22", "23"}, + ] + actual = [set(restriction) for restriction in mapk_cascade.restrictions] + for l1 in actual: + assert l1 in desired + for l2 in desired: + assert l2 in actual + + with pytest.raises(DuplicateError): + Text2Model( + os.path.join(os.path.dirname(__file__), "text_files", "duplicate_binding.txt") + ).convert() + + +def test_run_simulation(): + try: + from .text_files import Kholodenko1999, michaelis_menten + + _packaging: Callable[[str], str] = lambda model_name: ".".join( + ["tests.test_text2model.text_files", model_name] + ) + for model_name in ["Kholodenko1999", "michaelis_menten"]: + model = Model(_packaging(model_name)).create() + run_simulation(model) + simulated_values = np.load( + os.path.join( + os.path.dirname(__file__), + "text_files", + model_name, + "simulation_data", + "simulations_original.npy", + ) + ) + assert np.isfinite(simulated_values).all() + if model_name == "Kholodenko1999": + expected = np.load( + os.path.join( + os.path.dirname(__file__), + "simulations_original_BN.npy", + ) + ) + assert np.allclose(simulated_values, expected) + except ImportError as e: + print(e) + + +def test_text2markdown(): + for model in ["michaelis_menten", "Kholodenko1999"]: + if model == "michaelis_menten": + mm_kinetics = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"{model}.txt", + ) + ) + mm_kinetics.register_word({"dissociate": ["releases"]}) + mm_kinetics.to_markdown(num_reactions=2) + elif model == "Kholodenko1999": + mapk_cascade = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"{model}.txt", + ) + ) + mapk_cascade.to_markdown(num_reactions=25) + assert os.path.isfile(os.path.join("markdown", model, "rate.md")) + assert os.path.isfile(os.path.join("markdown", model, "diffeq.md")) + shutil.rmtree("markdown") + + +def test_julia_models(): + necessities = [ + os.path.join("name2idx", "parameters.jl"), + os.path.join("name2idx", "species.jl"), + "ode.jl", + "observable.jl", + "simulation.jl", + "experimental_data.jl", + "search_param.jl", + "problem.jl", + ] + for model in ["michaelis_menten", "Kholodenko1999"]: + for file in necessities: + assert os.path.isfile( + os.path.join( + os.path.dirname(__file__), + "text_files", + f"{model}_jl", + file, + ) + ) + + +def test_cleanup(): + for model in [ + "michaelis_menten", + "Kholodenko1999", + "michaelis_menten_jl", + "Kholodenko1999_jl", + "duplicate_binding", + ]: + assert os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + model, + ) + ) + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + model, + ) + ) diff --git a/tests/test_text2model/test_thermodynamic_restriction.py b/tests/test_text2model/test_thermodynamic_restriction.py new file mode 100644 index 00000000..08ff5f03 --- /dev/null +++ b/tests/test_text2model/test_thermodynamic_restriction.py @@ -0,0 +1,54 @@ +import os +import shutil + +from biomass import Text2Model + + +def test_preprocessing(): + if os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + "abc", + ) + ): + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + "abc", + ) + ) + + +def test_text2model(): + abc = Text2Model( + os.path.join( + os.path.dirname(__file__), + "text_files", + "abc.txt", + ), + ) + abc.convert() + # test thermodynamic restrictions + desired = [["1", "2", "3", "4"]] + actual = abc.restrictions + assert len(actual) == 1 + assert set(desired[0]) == set(actual[0]) + + +def test_cleanup(): + assert os.path.isdir( + os.path.join( + os.path.dirname(__file__), + "text_files", + "abc", + ) + ) + shutil.rmtree( + os.path.join( + os.path.dirname(__file__), + "text_files", + "abc", + ) + ) diff --git a/tests/test_text2model/text_files/Kholodenko1999.txt b/tests/test_text2model/text_files/Kholodenko1999.txt new file mode 100644 index 00000000..b09eb227 --- /dev/null +++ b/tests/test_text2model/text_files/Kholodenko1999.txt @@ -0,0 +1,44 @@ +EGF binds EGFR <--> Ra | kf=0.003, kr=0.06 | EGFR=100 +Ra dimerizes <--> R2 | kf=0.01, kr=0.1 +R2 is phosphorylated <--> RP | kf=1, kr=0.01 +RP is dephosphorylated --> R2 | V=450, K=50 +RP binds PLCg <--> RPL | kf=0.06, kr=0.2 | PLCg=105 +RPL is phosphorylated <--> RPLP | kf=1, kr=0.05 +RPLP is dissociated into RP and PLCgP | kf=0.3, kr=0.006 +PLCgP is dephosphorylated --> PLCg | V=1, K=100 +RP binds Grb2 <--> RG | kf=0.003, kr=0.05 | Grb2=85 +RG binds SOS <--> RGS | kf=0.01, kr=0.06 | SOS=34 +RGS is dissociated into RP and GS | kf=0.03, kr=4.5e-3 +GS is dissociated into Grb2 and SOS | kf=1.5e-3, kr=1e-4 +RP binds Shc <--> RSh | kf=0.09, kr=0.6 | Shc=150 +RSh is phosphorylated <--> RShP | kf=6, kr=0.06 +RShP is dissociated into ShP and RP | kf=0.3, kr=9e-4 +ShP is dephosphorylated --> Shc | V=1.7, K=340 +RShP binds Grb2 <--> RShG | kf=0.003, kr=0.1 +RShG is dissociated into RP and ShG | kf=0.3, kr=9e-4 +RShG binds SOS <--> RShGS | kf=0.01, kr=2.14e-2 +RShGS is dissociated into ShGS and RP | kf=0.12, kr=2.4e-4 +ShP binds Grb2 <--> ShG | kf=0.003, kr=0.1 +ShG binds SOS <--> ShGS | kf=0.03, kr=0.064 +ShGS is dissociated into ShP and GS | kf=0.1, kr=0.021 +RShP binds GS <--> RShGS | kf=0.009, kr=4.29e-2 +PLCgP is translocated to cytoskeletal or membrane structures <--> PLCgP_I | kf=1, kr=0.03 + +# Kholodenko, B. N., Demin, O. V, Moehren, G. & Hoek, J. B. +# Quantification of short term signaling by the epidermal growth factor receptor. +# J. Biol. Chem. 274, 30169–30181 (1999). https://doi.org/10.1074/jbc.274.42.30169 + +# observable layer +@obs Total_phosphorylated_Shc: u[RShP] + u[RShG] + u[RShGS] + u[ShP] + u[ShG] + u[ShGS] +@obs Total_Grb2_coprecipitated_with_Shc: u[RShG] + u[ShG] + u[RShGS] + u[ShGS] +@obs Total_phosphorylated_Shc_bound_to_EGFR: u[RShP] + u[RShG] + u[RShGS] +@obs Total_Grb2_bound_to_EGFR: u[RG] + u[RGS] + u[RShG] + u[RShGS] +@obs Total_SOS_bound_to_EGFR: u[RGS] + u[RShGS] +@obs ShGS_complex: u[ShGS] +@obs Total_phosphorylated_PLCg: u[RPLP] + u[PLCgP] + +# simulation layer +@sim tspan: [0, 120] +@sim condition EGF20nM: init[EGF] = 680 +@sim condition EGF2nM: init[EGF] = 68 +@sim condition Absence_PLCgP_transloc: init[EGF] = 680; p[kf25] = 0; p[kr25] = 0 \ No newline at end of file diff --git a/tests/test_text2model/text_files/abc.txt b/tests/test_text2model/text_files/abc.txt new file mode 100644 index 00000000..5b6ca8c0 --- /dev/null +++ b/tests/test_text2model/text_files/abc.txt @@ -0,0 +1,4 @@ +A + B <--> AB +AB + C <--> ABC +ABC <--> A + BC +BC <--> B + C \ No newline at end of file diff --git a/tests/test_text2model/text_files/duplicate_binding.txt b/tests/test_text2model/text_files/duplicate_binding.txt new file mode 100644 index 00000000..d5a83974 --- /dev/null +++ b/tests/test_text2model/text_files/duplicate_binding.txt @@ -0,0 +1,2 @@ +A + B <--> AB +AB <--> A + B \ No newline at end of file diff --git a/tests/test_text2model/text_files/fos_model.txt b/tests/test_text2model/text_files/fos_model.txt new file mode 100644 index 00000000..493f90ec --- /dev/null +++ b/tests/test_text2model/text_files/fos_model.txt @@ -0,0 +1,80 @@ +@rxn ERKc --> pERKc: p[V1] * p[a] * u[ppMEKc] * u[ERKc] / ( p[K1] * (1 + u[pERKc] / p[K2]) + u[ERKc] ) || ERKc=9.60e02 +@rxn pERKc --> ppERKc: p[V2] * p[a] * u[ppMEKc] * u[pERKc] / ( p[K2] * (1 + u[ERKc] / p[K1]) + u[pERKc] ) | const V2=2.20e-01, const K2=3.50e02 +@rxn pERKc --> ERKc: p[V3] * u[pERKc] / ( p[K3] * (1 + u[ppERKc] / p[K4]) + u[pERKc] ) | const V3=7.20e-01, const K3=1.60e02 +@rxn ppERKc --> pERKc: p[V4] * u[ppERKc] / ( p[K4]* (1 + u[pERKc] / p[K3]) + u[ppERKc] ) | const V4=6.48e-01, const K4=6.00e01 +@rxn pERKn --> ERKn: p[V5] * u[pERKn] / ( p[K5] * (1 + u[ppERKn] / p[K6]) + u[pERKn] ) +@rxn ppERKn --> pERKn: p[V6] * u[ppERKn] / ( p[K6] * (1 + u[pERKn] / p[K5]) + u[ppERKn] ) |5| +ERKc translocates to nucleus (0.94, 0.22) <--> ERKn | const kf=1.20e-02, const kr=1.80e-02 +pERKc translocates to nucleus (0.94, 0.22) <--> pERKn | const kf=1.20e-02, const kr=1.80e-02 +ppERKc translocates to nucleus (0.94, 0.22) <--> ppERKn | const kf=1.10e-02, const kr=1.30e-02 +ppERKn transcribes PreduspmRNAn +PreduspmRNAn translocates to cytoplasm --> duspmRNAc +duspmRNAc is degraded +duspmRNAc is translated into DUSPc +ppERKc phosphorylates DUSPc --> pDUSPc +pDUSPc is dephosphorylated --> DUSPc +DUSPc is degraded | const kf=2.57e-04 +pDUSPc is degraded | const kf=9.63e-05 +DUSPc translocates to nucleus (0.94, 0.22) <--> DUSPn +pDUSPc translocates to nucleus (0.94, 0.22) <--> pDUSPn |18| +ppERKn phosphorylates DUSPn --> pDUSPn +pDUSPn is dephosphorylated --> DUSPn +DUSPn is degraded | const kf=2.57e-04 +pDUSPn is degraded | const kf=9.63e-05 +ppERKc phosphorylates RSKc --> pRSKc || RSKc=3.53e02 +pRSKc is dephosphorylated --> RSKc +pRSKc translocates to nucleus (0.94, 0.22) <--> pRSKn +pRSKn phosphorylates CREBn --> pCREBn || CREBn=1.00e03 +pCREBn is dephosphorylated --> CREBn +ppERKn phosphorylates Elk1n --> pElk1n || Elk1n=1.51e03 +pElk1n is dephosphorylated --> Elk1n +pCREBn & pElk1n transcribes PrecfosmRNAn, repressed by Fn +PrecfosmRNAn translocates to cytoplasm --> cfosmRNAc +cfosmRNAc is degraded +cfosmRNAc is translated into cFOSc +ppERKc phosphorylates cFOSc --> pcFOSc +pRSKc phosphorylates cFOSc --> pcFOSc +pcFOSc is dephosphorylated --> cFOSc +cFOSc is degraded | const kf=2.57e-04 +pcFOSc is degraded | const kf=9.63e-05 +cFOSc translocates to nucleus (0.94, 0.22) <--> cFOSn +pcFOSc translocates to nucleus (0.94, 0.22) <--> pcFOSn |40| +ppERKn phosphorylates cFOSn --> pcFOSn +pRSKn phosphorylates cFOSn --> pcFOSn +pcFOSn is dephosphorylated --> cFOSn +cFOSn is degraded | const kf=2.57e-04 +pcFOSn is degraded | const kf=9.63e-05 +DUSPn + ppERKn <--> DUSPn_ppERKn +DUSPn_ppERKn --> DUSPn + pERKn +DUSPn + pERKn <--> DUSPn_pERKn +DUSPn_pERKn --> DUSPn + ERKn +DUSPn + ERKn <--> DUSPn_ERKn +pDUSPn + ppERKn <--> pDUSPn_ppERKn |47| +pDUSPn_ppERKn --> pDUSPn + pERKn |48| +pDUSPn + pERKn <--> pDUSPn_pERKn |49| +pDUSPn_pERKn --> pDUSPn + ERKn |50| +pDUSPn + ERKn <--> pDUSPn_ERKn |51| +pcFOSn transcribes PreFmRNAn +PreFmRNAn translocates to cytoplasm --> FmRNAc +FmRNAc is degraded +FmRNAc is translated into Fc +Fc is degraded +Fc translocates to nucleus (0.94, 0.22) <--> Fn +Fn is degraded + +@add species ppMEKc +@add param Ligand + +@obs Phosphorylated_MEKc: u[ppMEKc] +@obs Phosphorylated_ERKc: u[pERKc] + u[ppERKc] +@obs Phosphorylated_RSKw: u[pRSKc] + u[pRSKn] * (0.22 / 0.94) +@obs Phosphorylated_CREBw: u[pCREBn] * (0.22 / 0.94) +@obs dusp_mRNA: u[duspmRNAc] +@obs cfos_mRNA: u[cfosmRNAc] +@obs cFos_Protein: (u[pcFOSn] + u[cFOSn]) * (0.22 / 0.94) + u[cFOSc] + u[pcFOSc] +@obs Phosphorylated_cFos: u[pcFOSn] * (0.22 / 0.94) + u[pcFOSc] + +@sim tspan: [0, 5400] +@sim unperturbed: p[Ligand] = 0 +@sim condition EGF: p[Ligand] = 1 +@sim condition HRG: p[Ligand] = 2 \ No newline at end of file diff --git a/tests/test_text2model/text_files/michaelis_menten.txt b/tests/test_text2model/text_files/michaelis_menten.txt new file mode 100644 index 00000000..af51bc31 --- /dev/null +++ b/tests/test_text2model/text_files/michaelis_menten.txt @@ -0,0 +1,10 @@ +E binds S <--> ES | kf=0.003, kr=0.001 | E=100, S=50 +ES releases E and P | kf=0.002, kr=0 + +@obs Substrate: u[S] +@obs E_free: u[E] +@obs E_total: u[E] + u[ES] +@obs Product: u[P] +@obs Complex: u[ES] + +@sim tspan: [0, 100] \ No newline at end of file diff --git a/tests/test_text2model/text_files/typo_1.txt b/tests/test_text2model/text_files/typo_1.txt new file mode 100644 index 00000000..93cc435f --- /dev/null +++ b/tests/test_text2model/text_files/typo_1.txt @@ -0,0 +1 @@ +A bindss B <--> AB \ No newline at end of file diff --git a/tests/test_text2model/text_files/typo_2.txt b/tests/test_text2model/text_files/typo_2.txt new file mode 100644 index 00000000..be596aaa --- /dev/null +++ b/tests/test_text2model/text_files/typo_2.txt @@ -0,0 +1 @@ +A bnds B <--> AB \ No newline at end of file