Repository: pasmopy/pasmopy (language: python)

Issue:
Proposal: v0.4 is compatible with `biomass` version 0.8.0 or more
- [x] Drop Python 3.7 support
- [x] Remove `Text2Model` library as it is supported in biomass
- [x] Update docs, tests


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 6d48ab7..8a28721 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -10,7 +10,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.7', '3.10']
+        python-version: ['3.8', '3.10']
 
     steps:
       - uses: actions/checkout@v2
@@ -18,6 +18,10 @@ jobs:
         uses: actions/setup-python@v2
         with:
           python-version: ${{ matrix.python-version }}
+      - name: Update apt-get
+        run: sudo apt-get update
+      - name: Install graphviz
+        run: sudo apt-get install graphviz graphviz-dev
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
diff --git a/README.md b/README.md
index 1db50d5..4720d04 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,9 @@
-# Pasmopy – Patient-Specific Modeling in Python
-
-[![overview](https://raw.githubusercontent.com/pasmopy/pasmopy/master/docs/_static/img/overview.png)](https://pasmopy.readthedocs.io)
+<br>
+<p align="center">
+    <a href="https://pasmopy.readthedocs.io/en/latest">
+        <img src="https://raw.githubusercontent.com/pasmopy/pasmopy/master/docs/_static/img/pasmopy-project-logo.png" width="400">
+    </a>
+</p>
 
 [![PyPI version](https://img.shields.io/pypi/v/pasmopy.svg?logo=PyPI&logoColor=white)](https://pypi.python.org/pypi/pasmopy)
 [![Actions Status](https://github.com/pasmopy/pasmopy/workflows/Tests/badge.svg)](https://github.com/pasmopy/pasmopy/actions)
@@ -31,6 +34,7 @@ You can also build the documentation locally by running the following commands:
 $ cd docs
 $ make html
 ```
+
 The site will live in `_build/html/index.html`.
 
 ## Installation
@@ -41,7 +45,9 @@ The latest stable release (and required dependencies) can be installed from [PyP
 $ pip install pasmopy
 ```
 
-Pasmopy requires Python 3.7+ to run.
+Pasmopy requires Python 3.8+ to run.
+
+:ledger: **Note: You will need to manually install [Graphviz](https://www.graphviz.org) (version 2.42 or later).**
 
 ## Example
 
@@ -100,6 +106,7 @@ Using Pasmopy, we built a mechanistic model of ErbB receptor signaling network,
 ## Reference
 
 - Imoto, H., Yamashiro, S. & Okada, M. A text-based computational framework for patient -specific modeling for classification of cancers. _iScience_ **25**, 103944 (2022). https://doi.org/10.1016/j.isci.2022.103944
+
 ## Author
 
 [Hiroaki Imoto](https://github.com/himoto)
diff --git a/docs/installation.rst b/docs/installation.rst
index 74d16c2..ffff057 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -1,7 +1,7 @@
 Installation
 ============
 
-Pasmopy requires Python 3.7+ to run.
+Pasmopy requires Python 3.8+ to run.
 
 Installing Pasmopy
 ------------------
diff --git a/docs/modules/construction.rst b/docs/modules/construction.rst
index 3af731a..46dd3d5 100644
--- a/docs/modules/construction.rst
+++ b/docs/modules/construction.rst
@@ -2,9 +2,6 @@ Construction of executable models from text (:py:mod:`pasmopy.construction`)
 ============================================================================
 
 .. toctree::
-   :maxdepth: 3
+   :maxdepth: 1
    
-   thermodynamic_restrictions
-   reaction_rules
-   text2model
-   template
\ No newline at end of file
+   text2model
\ No newline at end of file
diff --git a/docs/modules/reaction_rules.rst b/docs/modules/reaction_rules.rst
deleted file mode 100644
index 32a1bfb..0000000
--- a/docs/modules/reaction_rules.rst
+++ /dev/null
@@ -1,7 +0,0 @@
-Available reaction rules (:py:mod:`pasmopy.construction.reaction_rules`)
-------------------------------------------------------------------------
-
-.. autoclass:: pasmopy.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/modules/template.rst b/docs/modules/template.rst
deleted file mode 100644
index b4c7729..0000000
--- a/docs/modules/template.rst
+++ /dev/null
@@ -1,20 +0,0 @@
-Template for BioMASS model construction (:py:mod:`pasmopy.construction.template`)
----------------------------------------------------------------------------------
-
-.. automodule:: pasmopy.construction.template.ode
-   :members:
-
-.. automodule:: pasmopy.construction.template.observable
-   :members:
-   
-.. automodule:: pasmopy.construction.template.search_param
-   :members:
-
-.. automodule:: pasmopy.construction.template.problem
-   :members:
-
-.. automodule:: pasmopy.construction.template.viz
-   :members:
-
-.. automodule:: pasmopy.construction.template.reaction_network
-   :members:
\ No newline at end of file
diff --git a/docs/modules/text2model.rst b/docs/modules/text2model.rst
index 27dcee1..e6ba968 100644
--- a/docs/modules/text2model.rst
+++ b/docs/modules/text2model.rst
@@ -1,6 +1,2 @@
 Text-to-model conversion (:py:mod:`pasmopy.construction.text2model`)
 --------------------------------------------------------------------
-
-.. autoclass:: pasmopy.construction.text2model.Text2Model
-   :members: convert, to_markdown, register_word
-
diff --git a/docs/modules/thermodynamic_restrictions.rst b/docs/modules/thermodynamic_restrictions.rst
deleted file mode 100644
index 4a442d6..0000000
--- a/docs/modules/thermodynamic_restrictions.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-Thermodynamic restrictions (:py:mod:`pasmopy.construction.thermodynamic_restrictions`)
---------------------------------------------------------------------------------------
-
-.. autoclass:: pasmopy.construction.thermodynamic_restrictions.ThermodynamicRestrictions
-   :members:
-
diff --git a/pasmopy/construction/__init__.py b/pasmopy/construction/__init__.py
index 0c10b9c..88880ae 100644
--- a/pasmopy/construction/__init__.py
+++ b/pasmopy/construction/__init__.py
@@ -1,1 +1,1 @@
-from .text2model import Text2Model
+from .text2model import Text2Model, run_example_model
diff --git a/pasmopy/construction/julia_template.py b/pasmopy/construction/julia_template.py
deleted file mode 100644
index 0780eb6..0000000
--- a/pasmopy/construction/julia_template.py
+++ /dev/null
@@ -1,525 +0,0 @@
-"""
-Model template for conversion of biomass models into
-BioMASS.jl (https://github.com/biomass-dev/BioMASS.jl) format.
-
-BioMASS.jl >= 0.7.2
-"""
-
-PARAMETERS: 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: 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: 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: 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: 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: 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: 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: 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/pasmopy/construction/reaction_rules.py b/pasmopy/construction/reaction_rules.py
deleted file mode 100644
index abd03b3..0000000
--- a/pasmopy/construction/reaction_rules.py
+++ /dev/null
@@ -1,1596 +0,0 @@
-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 <https://github.com/pasmopy/pasmopy/blob/master/tests/text_files/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 <observable name>: <expression>`."
-                )
-            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/pasmopy/construction/template/README.md b/pasmopy/construction/template/README.md
deleted file mode 100644
index 732c122..0000000
--- a/pasmopy/construction/template/README.md
+++ /dev/null
@@ -1,15 +0,0 @@
-# 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/pasmopy/construction/template/__init__.py b/pasmopy/construction/template/__init__.py
deleted file mode 100644
index a0bc163..0000000
--- a/pasmopy/construction/template/__init__.py
+++ /dev/null
@@ -1,5 +0,0 @@
-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/pasmopy/construction/template/name2idx/__init__.py b/pasmopy/construction/template/name2idx/__init__.py
deleted file mode 100644
index 3869ddc..0000000
--- a/pasmopy/construction/template/name2idx/__init__.py
+++ /dev/null
@@ -1,2 +0,0 @@
-from .parameters import C
-from .species import V
diff --git a/pasmopy/construction/template/name2idx/parameters.py b/pasmopy/construction/template/name2idx/parameters.py
deleted file mode 100644
index cb3917c..0000000
--- a/pasmopy/construction/template/name2idx/parameters.py
+++ /dev/null
@@ -1,19 +0,0 @@
-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/pasmopy/construction/template/name2idx/species.py b/pasmopy/construction/template/name2idx/species.py
deleted file mode 100644
index 8109670..0000000
--- a/pasmopy/construction/template/name2idx/species.py
+++ /dev/null
@@ -1,19 +0,0 @@
-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/pasmopy/construction/template/observable.py b/pasmopy/construction/template/observable.py
deleted file mode 100644
index 5b1b796..0000000
--- a/pasmopy/construction/template/observable.py
+++ /dev/null
@@ -1,77 +0,0 @@
-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/pasmopy/construction/template/ode.py b/pasmopy/construction/template/ode.py
deleted file mode 100644
index 3475633..0000000
--- a/pasmopy/construction/template/ode.py
+++ /dev/null
@@ -1,34 +0,0 @@
-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/pasmopy/construction/template/problem.py b/pasmopy/construction/template/problem.py
deleted file mode 100644
index 90ef863..0000000
--- a/pasmopy/construction/template/problem.py
+++ /dev/null
@@ -1,103 +0,0 @@
-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/pasmopy/construction/template/reaction_network.py b/pasmopy/construction/template/reaction_network.py
deleted file mode 100644
index f5523ed..0000000
--- a/pasmopy/construction/template/reaction_network.py
+++ /dev/null
@@ -1,39 +0,0 @@
-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/pasmopy/construction/template/search_param.py b/pasmopy/construction/template/search_param.py
deleted file mode 100644
index 70124b8..0000000
--- a/pasmopy/construction/template/search_param.py
+++ /dev/null
@@ -1,77 +0,0 @@
-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/pasmopy/construction/template/viz.py b/pasmopy/construction/template/viz.py
deleted file mode 100644
index 54114e9..0000000
--- a/pasmopy/construction/template/viz.py
+++ /dev/null
@@ -1,78 +0,0 @@
-from typing import List
-
-from biomass.plotting import *
-from matplotlib import pyplot as plt
-
-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/pasmopy/construction/text2model.py b/pasmopy/construction/text2model.py
index d2dac01..65faee2 100644
--- a/pasmopy/construction/text2model.py
+++ b/pasmopy/construction/text2model.py
@@ -1,911 +1,42 @@
-import os
-import re
-import shutil
-from dataclasses import dataclass, field
-from typing import Dict, List, Optional
+"""
+:py:class:`Text2Model` is a useful class to build an ordinary differential equation (ODE) model from a text file describing biochemical systems.
+It was originally developed in the pasmopy project [1]_ but migrated to the `biomass <https://github.com/biomass-dev/biomass>`_ framework.
 
-try:  # python 3.8+
-    from typing import Literal
-except ImportError:
-    from typing_extensions import Literal
+* **How to use:** https://pasmopy.readthedocs.io/en/latest/model_development.html
+* **API reference:** https://biomass-core.readthedocs.io/en/latest/api/construction.html
 
-from . import julia_template as jl
-from .reaction_rules import ReactionRules
+References
+----------
+.. [1] Imoto, H., Yamashiro, S. & Okada, M.
+    A text-based computational framework for patient -specific modeling for classification of cancers.
+    iScience 25, 103944 (2022).
+"""
+from typing import Final
 
+from biomass import Text2Model, create_model, run_simulation
 
-@dataclass
-class Text2Model(ReactionRules):
-    """
-    Build a BioMASS-formatted model based on template.
-
-    **reaction** | **parameters** | **initial conditions**
+MICHAELIS_MENTEN: Final = """\n
+E + S ⇄ ES | kf=0.003, kr=0.001 | E=100, S=50
+ES → E + P | kf=0.002
 
-    Attributes
-    ----------
-    input_txt : str
-        Model description file (*.txt), e.g., 'Kholodenko1999.txt'
+@obs Substrate: u[S]
+@obs E_free: u[E]
+@obs E_total: u[E] + u[ES]
+@obs Product: u[P]
+@obs Complex: u[ES]
 
-    similarity_threshold : float (default: 0.7)
-        Similarity threshold used in text-to-model conversion. Must lie within (0, 1).
+@sim tspan: [0, 100]
+"""
 
-    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)
+def run_example_model() -> None:
     """
-
-    input_txt: str
-    similarity_threshold: float = 0.7
-    lang: Literal["python", "julia"] = "python"
-    indentation: 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, n_reaction: int, savedir: str = "markdown") -> None:
-        """
-        Create markdown table describing differential equations.
-
-        Parameters
-        ----------
-        n_reaction : 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 n_reaction <= 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"<sup>{p_names.strip()}</sup>".replace(
-                                f"{num + 1:d}", f"<sub>{num + 1:d}</sub>"
-                            ),
-                        )
-                    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"<sub>{num + 1:d}</sub>"
-                            ),
-                        )
-                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[:n_reaction],
-                ]
-            )
-        # 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(n_reaction):
-                if f"v[{num + 1:d}]" in eq:
-                    eq = eq.replace(
-                        f"v[{num + 1:d}]",
-                        f"_v_ <sub>{num + 1:d}</sub>",
-                    )
-            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)
+    A simple Michaelis-Menten two-step enzyme catalysis model.
+    See README at https://github.com/pasmopy/pasmopy.
+    """
+    with open("michaelis_menten.txt", mode="w") as f:
+        f.writelines(MICHAELIS_MENTEN)
+    description = Text2Model("michaelis_menten.txt")
+    description.convert()
+    model = create_model("michaelis_menten")
+    run_simulation(model)
diff --git a/pasmopy/construction/thermodynamic_restrictions.py b/pasmopy/construction/thermodynamic_restrictions.py
deleted file mode 100644
index 3ce5363..0000000
--- a/pasmopy/construction/thermodynamic_restrictions.py
+++ /dev/null
@@ -1,205 +0,0 @@
-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/pasmopy/patient_model.py b/pasmopy/patient_model.py
index fbd7b7a..77b6276 100644
--- a/pasmopy/patient_model.py
+++ b/pasmopy/patient_model.py
@@ -2,12 +2,7 @@
 import multiprocessing
 import os
 from dataclasses import dataclass, field
-from typing import Callable, Dict, List, Optional, Union
-
-try:  # python 3.8+
-    from typing import Literal
-except ImportError:
-    from typing_extensions import Literal
+from typing import Callable, Dict, List, Literal, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -82,8 +77,7 @@ def _check_ctx(context: str) -> None:
         """
         Check whether the context is appropriate.
         """
-        contexts = ["spawn", "fork", "forkserver"]
-        if context not in contexts:
+        if context not in (contexts := ["spawn", "fork", "forkserver"]):
             raise ValueError("context must be one of '{}'.".format("', '".join(contexts)))
 
 
diff --git a/pasmopy/version.py b/pasmopy/version.py
index 260c070..6a9beea 100644
--- a/pasmopy/version.py
+++ b/pasmopy/version.py
@@ -1,1 +1,1 @@
-__version__ = "0.3.1"
+__version__ = "0.4.0"
diff --git a/pyproject.toml b/pyproject.toml
index fc4eb8a..18e3ae7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.black]
 line-length = 99
-target-version = ['py37']
+target-version = ['py39']
 exclude = '''
 /(
     \.eggs
diff --git a/requirements.txt b/requirements.txt
index 81c96ee..41a905b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,6 @@
-biomass>=0.7
+biomass>=0.8
 numpy>=1.17
 pandas>=0.24
 seaborn>=0.11.2
 scipy>=1.6
-tqdm>=4.50.2
-typing_extensions
\ No newline at end of file
+tqdm>=4.50.2
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 2f03493..bf97c35 100644
--- a/setup.py
+++ b/setup.py
@@ -22,8 +22,8 @@ def get_version() -> str:
 
 def setup_package():
     # Python version check.
-    if sys.version_info[:2] < (3, 7):
-        raise RuntimeError("Pasmopy requires at least Python version 3.7")
+    if sys.version_info[:2] < (3, 8):
+        raise RuntimeError("Pasmopy requires at least Python version 3.8")
 
     setup(
         name="pasmopy",
@@ -62,7 +62,7 @@ def setup_package():
                 "sphinxcontrib-bibtex>=2.2",
             ],
         },
-        python_requires=">=3.7",
+        python_requires=">=3.8",
         keywords=[
             "systems",
             "biology",
@@ -82,7 +82,6 @@ def setup_package():
             "Natural Language :: English",
             "Operating System :: OS Independent",
             "Programming Language :: Python :: 3",
-            "Programming Language :: Python :: 3.7",
             "Programming Language :: Python :: 3.8",
             "Programming Language :: Python :: 3.9",
             "Programming Language :: Python :: 3.10",


Test patch:
diff --git a/tests/C.py b/tests/C.py
index 5306da9..216cbfd 100755
--- a/tests/C.py
+++ b/tests/C.py
@@ -1,240 +1,3 @@
-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
-"""
-
-
 REQUIREMENTS = """\
 import os
 from pasmopy import Individualization
@@ -277,12 +40,12 @@ def get_timepoint(obs_name) -> List[int]:
 
 INCORPORATION = """\
         x[C.V291] = incorporating_gene_expression_levels.as_reaction_rate(
-            __path__[0].split(os.sep)[-1].replace("_", "."), x, "V291", "DUSP"
+            __path__[0].split(os.sep)[-1], x, "V291", "DUSP"
         )
         x[C.V310] = incorporating_gene_expression_levels.as_reaction_rate(
-            __path__[0].split(os.sep)[-1].replace("_", "."), x, "V310", "cMyc"
+            __path__[0].split(os.sep)[-1], x, "V310", "cMyc"
         )
         y0 = incorporating_gene_expression_levels.as_initial_conditions(
-            __path__[0].split(os.sep)[-1].replace("_", "."), x, y0
+            __path__[0].split(os.sep)[-1], x, y0
         )
 """
diff --git a/tests/models/breast/sample_names.txt b/tests/models/breast/sample_names.txt
index 5d5c658..1dc3d10 100644
--- a/tests/models/breast/sample_names.txt
+++ b/tests/models/breast/sample_names.txt
@@ -1,377 +1,369 @@
-TCGA_3C_AALK_01A
-TCGA_A1_A0SD_01A
-TCGA_A1_A0SE_01A
-TCGA_A1_A0SF_01A
-TCGA_A1_A0SH_01A
-TCGA_A1_A0SI_01A
-TCGA_A1_A0SK_01A
-TCGA_A1_A0SN_01A
+TCGA_A2_A0D4_01A
+TCGA_D8_A143_01A
+TCGA_C8_A134_01A
+TCGA_E9_A1NA_01A
+TCGA_BH_A28Q_01A
+TCGA_E2_A1IO_01A
+TCGA_GM_A2DF_01A
+TCGA_A2_A25B_01A
+TCGA_BH_A0B0_01A
+TCGA_E2_A14V_01A
+TCGA_AN_A0FF_01A
+TCGA_AC_A2BM_01A
 TCGA_A1_A0SP_01A
-TCGA_A1_A0SQ_01A
-TCGA_A2_A04Q_01A
-TCGA_A2_A04R_01A
-TCGA_A2_A04U_01A
-TCGA_A2_A04V_01A
-TCGA_A2_A04W_01A
-TCGA_A2_A04X_01A
-TCGA_A2_A04Y_01A
-TCGA_A2_A0CM_01A
+TCGA_BH_A6R8_01A
+TCGA_A8_A0AB_01A
+TCGA_A7_A26G_01A
+TCGA_BH_A0B7_01A
+TCGA_E9_A244_01A
+TCGA_EW_A1P4_01A
+TCGA_D8_A1JD_01A
+TCGA_BH_A18K_01A
+TCGA_D8_A27F_01A
+TCGA_A7_A0D9_01A
+TCGA_E2_A15K_01A
+TCGA_E2_A15K_06A
+TCGA_E2_A14Y_01A
+TCGA_B6_A0I5_01A
+TCGA_A2_A3XV_01A
+TCGA_A2_A0SX_01A
+TCGA_AR_A24Z_01A
+TCGA_E2_A15H_01A
+TCGA_AC_A7VB_01A
+TCGA_AO_A0J2_01A
+TCGA_AR_A256_01A
+TCGA_OL_A6VR_01A
+TCGA_BH_A0DG_01A
+TCGA_E2_A1BD_01A
+TCGA_AO_A03T_01A
+TCGA_BH_A2L8_01A
+TCGA_D8_A1JI_01A
+TCGA_AR_A1AH_01A
+TCGA_E2_A153_01A
+TCGA_AR_A1AN_01A
+TCGA_OL_A5RZ_01A
+TCGA_GM_A2DN_01A
+TCGA_BH_A42V_01A
+TCGA_D8_A1XU_01A
+TCGA_E9_A22G_01A
+TCGA_E2_A1IG_01A
+TCGA_AR_A2LM_01A
+TCGA_AR_A0TS_01A
+TCGA_BH_A0HA_01A
+TCGA_D8_A1XL_01A
+TCGA_AO_A0JJ_01A
 TCGA_A2_A0CR_01A
+TCGA_A8_A095_01A
+TCGA_BH_A1FE_01A
+TCGA_A2_A0EO_01A
+TCGA_A2_A3XT_01A
+TCGA_A2_A1G0_01A
+TCGA_A8_A0A7_01A
+TCGA_BH_A0BF_01A
+TCGA_D8_A146_01A
+TCGA_3C_AALK_01A
 TCGA_A2_A0CV_01A
-TCGA_A2_A0CX_01A
-TCGA_A2_A0CZ_01A
+TCGA_EW_A6S9_01A
+TCGA_BH_A0BL_01A
+TCGA_EW_A1J5_01A
+TCGA_EW_A1IX_01A
+TCGA_E2_A1L6_01A
+TCGA_E2_A1LB_01A
+TCGA_EW_A2FW_01A
+TCGA_E2_A1B5_01A
+TCGA_E2_A14N_01A
 TCGA_A2_A0D2_01A
-TCGA_A2_A0D3_01A
-TCGA_A2_A0D4_01A
-TCGA_A2_A0EO_01A
-TCGA_A2_A0EP_01A
-TCGA_A2_A0ES_01A
+TCGA_E2_A10C_01A
 TCGA_A2_A0EX_01A
-TCGA_A2_A0SX_01A
-TCGA_A2_A0T0_01A
-TCGA_A2_A0T3_01A
-TCGA_A2_A0T5_01A
-TCGA_A2_A0T6_01A
-TCGA_A2_A0T7_01A
-TCGA_A2_A0YC_01A
-TCGA_A2_A0YE_01A
-TCGA_A2_A1G0_01A
-TCGA_A2_A25A_01A
-TCGA_A2_A25B_01A
-TCGA_A2_A25C_01A
-TCGA_A2_A3KC_01A
-TCGA_A2_A3XT_01A
+TCGA_BH_A0DT_01A
+TCGA_A7_A0CE_01A
+TCGA_AR_A24W_01A
+TCGA_C8_A8HQ_01A
+TCGA_B6_A0RU_01A
+TCGA_BH_A18V_06A
+TCGA_HN_A2NL_01A
+TCGA_A2_A04W_01A
+TCGA_EW_A1P7_01A
+TCGA_BH_A0DQ_01A
+TCGA_C8_A1HJ_01A
 TCGA_A2_A3XU_01A
-TCGA_A2_A3XV_01A
-TCGA_A2_A3XW_01A
-TCGA_A2_A3XX_01A
+TCGA_A2_A0T3_01A
+TCGA_B6_A0IJ_01A
+TCGA_E2_A15E_06A
 TCGA_A2_A3XY_01A
-TCGA_A2_A3XZ_01A
+TCGA_A2_A3KC_01A
+TCGA_OL_A66O_01A
+TCGA_BH_A0BT_01A
+TCGA_E9_A248_01A
 TCGA_A2_A3Y0_01A
-TCGA_A7_A0CE_01A
+TCGA_LL_A6FR_01A
+TCGA_E9_A22H_01A
+TCGA_A1_A0SI_01A
+TCGA_AR_A250_01A
+TCGA_GI_A2C9_01A
+TCGA_A2_A0CZ_01A
+TCGA_E9_A243_01A
+TCGA_EW_A2FS_01A
 TCGA_A7_A0CJ_01A
-TCGA_A7_A0D9_01A
-TCGA_A7_A13D_01A
-TCGA_A7_A26G_01A
-TCGA_A7_A26J_01A
-TCGA_A7_A4SD_01A
-TCGA_A7_A5ZW_01A
-TCGA_A7_A6VV_01A
-TCGA_A7_A6VW_01A
-TCGA_A7_A6VY_01A
-TCGA_A8_A075_01A
-TCGA_A8_A07C_01A
-TCGA_A8_A07J_01A
-TCGA_A8_A07O_01A
-TCGA_A8_A082_01A
-TCGA_A8_A086_01A
-TCGA_A8_A08G_01A
+TCGA_BH_A0BJ_01A
+TCGA_E2_A159_01A
+TCGA_AO_A03V_01A
+TCGA_B6_A402_01A
+TCGA_E2_A10A_01A
 TCGA_A8_A08I_01A
+TCGA_GM_A2DM_01A
+TCGA_A2_A25C_01A
+TCGA_BH_A0C3_01A
+TCGA_A2_A0T7_01A
+TCGA_GM_A2DH_01A
+TCGA_D8_A1XF_01A
+TCGA_E9_A247_01A
+TCGA_EW_A6SD_01A
+TCGA_AR_A24Q_01A
+TCGA_AC_A2FF_01A
+TCGA_BH_A0EI_01A
 TCGA_A8_A08R_01A
-TCGA_A8_A095_01A
-TCGA_A8_A09D_01A
-TCGA_A8_A09V_01A
-TCGA_A8_A0A7_01A
-TCGA_A8_A0AB_01A
+TCGA_A8_A082_01A
+TCGA_E2_A1B1_01A
+TCGA_C8_A1HN_01A
 TCGA_AC_A23E_01A
-TCGA_AC_A2BM_01A
-TCGA_AC_A2FF_01A
+TCGA_AR_A1AS_01A
+TCGA_D8_A27K_01A
+TCGA_AR_A24K_01A
+TCGA_A1_A0SK_01A
+TCGA_E2_A574_01A
 TCGA_AC_A6IV_01A
-TCGA_AC_A6NO_01A
-TCGA_AC_A7VB_01A
-TCGA_AC_A7VC_01A
-TCGA_AN_A041_01A
+TCGA_GM_A2DK_01A
+TCGA_E9_A5UO_01A
+TCGA_BH_A1ES_06A
+TCGA_A2_A04Q_01A
+TCGA_A8_A08G_01A
+TCGA_AR_A0U4_01A
+TCGA_B6_A0RG_01A
+TCGA_BH_A18I_01A
 TCGA_AN_A04C_01A
-TCGA_AN_A04D_01A
-TCGA_AN_A0AM_01A
-TCGA_AN_A0AR_01A
-TCGA_AN_A0FF_01A
-TCGA_AN_A0FS_01A
-TCGA_AN_A0FV_01A
-TCGA_AN_A0FX_01A
-TCGA_AN_A0FY_01A
-TCGA_AN_A0G0_01A
+TCGA_GM_A3XL_01A
+TCGA_B6_A0I6_01A
+TCGA_C8_A1HF_01A
+TCGA_BH_A5IZ_01A
+TCGA_A2_A0ES_01A
+TCGA_D8_A1JM_01A
 TCGA_AN_A0XU_01A
-TCGA_AO_A03P_01A
-TCGA_AO_A03T_01A
-TCGA_AO_A03V_01A
-TCGA_AO_A0J2_01A
-TCGA_AO_A0J4_01A
-TCGA_AO_A0JI_01A
-TCGA_AO_A0JJ_01A
-TCGA_AO_A124_01A
-TCGA_AO_A126_01A
-TCGA_AO_A129_01A
-TCGA_AO_A12A_01A
-TCGA_AO_A12C_01A
-TCGA_AO_A12D_01A
-TCGA_AO_A12E_01A
-TCGA_AO_A12F_01A
-TCGA_AO_A1KO_01A
-TCGA_AO_A1KR_01A
-TCGA_AQ_A04J_01A
-TCGA_AR_A0TP_01A
-TCGA_AR_A0TS_01A
-TCGA_AR_A0TU_01A
-TCGA_AR_A0TY_01A
-TCGA_AR_A0U4_01A
-TCGA_AR_A1AH_01A
-TCGA_AR_A1AI_01A
-TCGA_AR_A1AN_01A
-TCGA_AR_A1AO_01A
-TCGA_AR_A1AQ_01A
-TCGA_AR_A1AS_01A
-TCGA_AR_A24K_01A
+TCGA_JL_A3YW_01A
+TCGA_AN_A0FY_01A
 TCGA_AR_A24L_01A
-TCGA_AR_A24N_01A
+TCGA_BH_A8FZ_01A
+TCGA_A7_A5ZW_01A
+TCGA_AN_A0FX_01A
+TCGA_BH_A0B3_01A
+TCGA_D8_A13Y_01A
+TCGA_AC_A6NO_01A
+TCGA_A2_A04V_01A
 TCGA_AR_A24P_01A
-TCGA_AR_A24Q_01A
-TCGA_AR_A24U_01A
-TCGA_AR_A24V_01A
-TCGA_AR_A24W_01A
-TCGA_AR_A24Z_01A
-TCGA_AR_A250_01A
-TCGA_AR_A252_01A
-TCGA_AR_A256_01A
-TCGA_AR_A2LM_01A
-TCGA_AR_A2LQ_01A
-TCGA_AR_A2LR_01A
-TCGA_AR_A5QP_01A
-TCGA_B6_A0I2_01A
-TCGA_B6_A0I5_01A
-TCGA_B6_A0I6_01A
-TCGA_B6_A0IA_01A
-TCGA_B6_A0IG_01A
-TCGA_B6_A0IJ_01A
-TCGA_B6_A0RG_01A
-TCGA_B6_A0RH_01A
-TCGA_B6_A0RS_01A
-TCGA_B6_A0RU_01A
-TCGA_B6_A0WZ_01A
-TCGA_B6_A401_01A
-TCGA_B6_A402_01A
-TCGA_B6_A40C_01A
-TCGA_BH_A0AU_01A
+TCGA_D8_A1J9_01A
+TCGA_BH_A18Q_01A
+TCGA_BH_A0HB_01A
+TCGA_LD_A7W6_01A
+TCGA_OL_A5D8_01A
+TCGA_BH_A1FN_01A
+TCGA_BH_A0BO_01A
+TCGA_D8_A1JB_01A
+TCGA_JL_A3YX_01A
 TCGA_BH_A0AV_01A
-TCGA_BH_A0AW_01A
-TCGA_BH_A0B0_01A
-TCGA_BH_A0B3_01A
-TCGA_BH_A0B6_01A
-TCGA_BH_A0B7_01A
+TCGA_A8_A075_01A
+TCGA_BH_A8G0_01A
 TCGA_BH_A0B9_01A
+TCGA_E2_A1LI_01A
+TCGA_OL_A5RY_01A
+TCGA_D8_A1XR_01A
+TCGA_LL_A442_01A
+TCGA_D8_A1JJ_01A
+TCGA_A2_A0EP_01A
+TCGA_E2_A155_01A
+TCGA_AR_A252_01A
+TCGA_A1_A0SN_01A
+TCGA_A7_A6VV_01A
+TCGA_BH_A1FE_06A
+TCGA_C8_A12Y_01A
+TCGA_A1_A0SH_01A
+TCGA_B6_A0IA_01A
+TCGA_A2_A3XX_01A
+TCGA_C8_A12V_01A
+TCGA_A2_A3XW_01A
+TCGA_BH_A0HL_01A
+TCGA_BH_A0W7_01A
 TCGA_BH_A0BD_01A
-TCGA_BH_A0BF_01A
-TCGA_BH_A0BJ_01A
-TCGA_BH_A0BL_01A
-TCGA_BH_A0BO_01A
+TCGA_EW_A1PH_01A
 TCGA_BH_A0BQ_01A
-TCGA_BH_A0BR_01A
-TCGA_BH_A0BT_01A
-TCGA_BH_A0C3_01A
-TCGA_BH_A0DG_01A
-TCGA_BH_A0DK_01A
-TCGA_BH_A0DQ_01A
-TCGA_BH_A0DT_01A
-TCGA_BH_A0E1_01A
-TCGA_BH_A0EI_01A
+TCGA_AO_A12A_01A
+TCGA_E9_A1N6_01A
+TCGA_E2_A15S_01A
+TCGA_A8_A09V_01A
+TCGA_E2_A15I_01A
+TCGA_C8_A26X_01A
+TCGA_C8_A12T_01A
+TCGA_A2_A04X_01A
+TCGA_AO_A12C_01A
+TCGA_BH_A1EW_01A
+TCGA_A2_A04Y_01A
+TCGA_D8_A1XO_01A
+TCGA_BH_A18V_01A
 TCGA_BH_A0H3_01A
-TCGA_BH_A0H5_01A
-TCGA_BH_A0HA_01A
-TCGA_BH_A0HB_01A
-TCGA_BH_A0HL_01A
-TCGA_BH_A0HO_01A
-TCGA_BH_A0HQ_01A
-TCGA_BH_A0HU_01A
-TCGA_BH_A0HX_01A
-TCGA_BH_A0RX_01A
-TCGA_BH_A0W3_01A
+TCGA_LL_A5YP_01A
+TCGA_AN_A0FS_01A
+TCGA_AN_A0AM_01A
+TCGA_OL_A5RW_01A
+TCGA_A2_A0YC_01A
 TCGA_BH_A0W4_01A
-TCGA_BH_A0W7_01A
-TCGA_BH_A18F_01A
-TCGA_BH_A18I_01A
-TCGA_BH_A18K_01A
-TCGA_BH_A18Q_01A
-TCGA_BH_A18R_01A
-TCGA_BH_A18V_06A
-TCGA_BH_A18V_01A
-TCGA_BH_A1ES_06A
-TCGA_BH_A1ET_01A
-TCGA_BH_A1EW_01A
-TCGA_BH_A1FE_01A
-TCGA_BH_A1FE_06A
-TCGA_BH_A1FN_01A
-TCGA_BH_A1FU_01A
-TCGA_BH_A208_01A
-TCGA_BH_A28Q_01A
-TCGA_BH_A2L8_01A
-TCGA_BH_A42V_01A
-TCGA_BH_A5IZ_01A
-TCGA_BH_A6R8_01A
-TCGA_BH_A8FZ_01A
-TCGA_BH_A8G0_01A
-TCGA_C8_A12N_01A
-TCGA_C8_A12O_01A
-TCGA_C8_A12P_01A
-TCGA_C8_A12T_01A
-TCGA_C8_A12U_01A
-TCGA_C8_A12V_01A
-TCGA_C8_A12Y_01A
-TCGA_C8_A12Z_01A
-TCGA_C8_A132_01A
-TCGA_C8_A134_01A
-TCGA_C8_A137_01A
-TCGA_C8_A1HE_01A
-TCGA_C8_A1HF_01A
-TCGA_C8_A1HJ_01A
-TCGA_C8_A1HK_01A
-TCGA_C8_A1HN_01A
-TCGA_C8_A26W_01A
-TCGA_C8_A26X_01A
+TCGA_A2_A0T6_01A
+TCGA_B6_A0WZ_01A
+TCGA_A1_A0SQ_01A
+TCGA_AO_A12E_01A
+TCGA_E2_A15J_01A
+TCGA_AC_A7VC_01A
 TCGA_C8_A26Z_01A
+TCGA_AR_A0TY_01A
+TCGA_BH_A0HU_01A
+TCGA_BH_A208_01A
+TCGA_OL_A5RX_01A
 TCGA_C8_A273_01A
-TCGA_C8_A275_01A
-TCGA_C8_A27A_01A
-TCGA_C8_A27B_01A
-TCGA_C8_A8HQ_01A
-TCGA_D8_A13Y_01A
-TCGA_D8_A143_01A
-TCGA_D8_A146_01A
-TCGA_D8_A1J9_01A
-TCGA_D8_A1JB_01A
-TCGA_D8_A1JD_01A
-TCGA_D8_A1JH_01A
-TCGA_D8_A1JI_01A
-TCGA_D8_A1JJ_01A
-TCGA_D8_A1JM_01A
-TCGA_D8_A1JU_01A
+TCGA_OK_A5Q2_01A
 TCGA_D8_A1X7_01A
-TCGA_D8_A1XF_01A
-TCGA_D8_A1XK_01A
-TCGA_D8_A1XL_01A
-TCGA_D8_A1XM_01A
-TCGA_D8_A1XO_01A
-TCGA_D8_A1XR_01A
-TCGA_D8_A1XU_01A
-TCGA_D8_A1XW_01A
-TCGA_D8_A27F_01A
-TCGA_D8_A27K_01A
-TCGA_D8_A27M_01A
-TCGA_D8_A73X_01A
-TCGA_E2_A106_01A
-TCGA_E2_A10A_01A
-TCGA_E2_A10C_01A
+TCGA_E2_A570_01A
+TCGA_A8_A07O_01A
+TCGA_E2_A1LA_01A
+TCGA_E9_A229_01A
+TCGA_E9_A22D_01A
 TCGA_E2_A10F_01A
-TCGA_E2_A14N_01A
+TCGA_E2_A1LS_01A
+TCGA_A2_A0T0_01A
+TCGA_E9_A1RB_01A
+TCGA_AN_A0G0_01A
+TCGA_OL_A6VO_01A
+TCGA_AR_A2LQ_01A
+TCGA_EW_A1PF_01A
+TCGA_AN_A0FV_01A
+TCGA_BH_A1ET_01A
 TCGA_E2_A14Q_01A
-TCGA_E2_A14T_01A
-TCGA_E2_A14V_01A
-TCGA_E2_A14Y_01A
-TCGA_E2_A150_01A
-TCGA_E2_A152_01A
-TCGA_E2_A153_01A
-TCGA_E2_A155_01A
-TCGA_E2_A158_01A
-TCGA_E2_A159_01A
-TCGA_E2_A15D_01A
-TCGA_E2_A15E_06A
-TCGA_E2_A15H_01A
-TCGA_E2_A15I_01A
-TCGA_E2_A15J_01A
-TCGA_E2_A15K_01A
-TCGA_E2_A15K_06A
-TCGA_E2_A15S_01A
-TCGA_E2_A1B1_01A
-TCGA_E2_A1B5_01A
-TCGA_E2_A1B6_01A
-TCGA_E2_A1BD_01A
-TCGA_E2_A1IG_01A
-TCGA_E2_A1IO_01A
-TCGA_E2_A1L6_01A
-TCGA_E2_A1L8_01A
+TCGA_AR_A24U_01A
+TCGA_AR_A24N_01A
+TCGA_BH_A0HX_01A
+TCGA_AO_A126_01A
+TCGA_AR_A24V_01A
+TCGA_Z7_A8R6_01A
+TCGA_D8_A1JU_01A
+TCGA_AO_A124_01A
+TCGA_BH_A0DK_01A
+TCGA_EW_A1J1_01A
+TCGA_D8_A1JH_01A
 TCGA_E2_A1L9_01A
-TCGA_E2_A1LA_01A
-TCGA_E2_A1LB_01A
-TCGA_E2_A1LG_01A
-TCGA_E2_A1LI_01A
-TCGA_E2_A1LS_01A
-TCGA_E2_A570_01A
-TCGA_E2_A573_01A
-TCGA_E2_A574_01A
-TCGA_E9_A1N5_01A
-TCGA_E9_A1N6_01A
-TCGA_E9_A1N8_01A
+TCGA_C8_A275_01A
+TCGA_AR_A2LR_01A
+TCGA_AR_A0TU_01A
+TCGA_EW_A1OZ_01A
+TCGA_E2_A1B6_01A
+TCGA_B6_A0RS_01A
+TCGA_D8_A1XK_01A
+TCGA_E9_A3X8_01A
+TCGA_OL_A66I_01A
+TCGA_HN_A2OB_01A
+TCGA_AR_A1AO_01A
+TCGA_BH_A0BR_01A
+TCGA_E9_A228_01A
 TCGA_E9_A1N9_01A
-TCGA_E9_A1NA_01A
-TCGA_E9_A1NE_01A
-TCGA_E9_A1NI_01A
-TCGA_E9_A1R0_01A
+TCGA_C8_A12U_01A
+TCGA_C8_A12N_01A
+TCGA_E9_A249_01A
+TCGA_B6_A40C_01A
+TCGA_AN_A041_01A
+TCGA_A2_A3XZ_01A
+TCGA_C8_A12O_01A
+TCGA_C8_A132_01A
+TCGA_B6_A0I2_01A
 TCGA_E9_A1RA_01A
-TCGA_E9_A1RB_01A
+TCGA_A1_A0SF_01A
 TCGA_E9_A227_01A
-TCGA_E9_A228_01A
-TCGA_E9_A229_01A
-TCGA_E9_A22D_01A
-TCGA_E9_A22G_01A
-TCGA_E9_A22H_01A
-TCGA_E9_A243_01A
-TCGA_E9_A244_01A
-TCGA_E9_A247_01A
-TCGA_E9_A248_01A
-TCGA_E9_A249_01A
-TCGA_E9_A3HO_01A
-TCGA_E9_A3QA_01A
-TCGA_E9_A3X8_01A
-TCGA_E9_A5UO_01A
-TCGA_EW_A1IX_01A
-TCGA_EW_A1IY_01A
-TCGA_EW_A1J1_01A
-TCGA_EW_A1J5_01A
-TCGA_EW_A1OV_01A
-TCGA_EW_A1OW_01A
+TCGA_PE_A5DE_01A
+TCGA_BH_A0RX_01A
+TCGA_A2_A04U_01A
+TCGA_AO_A0J4_01A
+TCGA_GM_A2DD_01A
+TCGA_AR_A0TP_01A
+TCGA_BH_A18R_01A
+TCGA_E2_A14T_01A
+TCGA_S3_AA14_01A
+TCGA_D8_A1XM_01A
+TCGA_BH_A0AU_01A
+TCGA_A8_A09D_01A
+TCGA_B6_A0IG_01A
+TCGA_A2_A0T5_01A
+TCGA_LD_A66U_01A
+TCGA_BH_A0E1_01A
 TCGA_EW_A1OX_01A
-TCGA_EW_A1OZ_01A
+TCGA_EW_A1IY_01A
+TCGA_E2_A573_01A
+TCGA_AR_A1AI_01A
+TCGA_C8_A12Z_01A
+TCGA_AR_A1AQ_01A
+TCGA_A7_A6VW_01A
 TCGA_EW_A1P0_01A
-TCGA_EW_A1P4_01A
-TCGA_EW_A1P7_01A
+TCGA_E9_A3HO_01A
+TCGA_AQ_A04J_01A
+TCGA_BH_A18F_01A
+TCGA_BH_A0HO_01A
+TCGA_A2_A0D3_01A
+TCGA_E2_A150_01A
+TCGA_A8_A07J_01A
+TCGA_AO_A12D_01A
+TCGA_D8_A73X_01A
+TCGA_C8_A12P_01A
+TCGA_A2_A25A_01A
+TCGA_AN_A0AR_01A
+TCGA_AO_A0JI_01A
+TCGA_E2_A106_01A
+TCGA_E9_A3QA_01A
+TCGA_E2_A1L8_01A
+TCGA_LL_A5YN_01A
+TCGA_C8_A26W_01A
+TCGA_A8_A086_01A
+TCGA_BH_A0B6_01A
+TCGA_S3_AA15_01A
+TCGA_A2_A0YE_01A
+TCGA_A7_A4SD_01A
+TCGA_AN_A04D_01A
 TCGA_EW_A1PA_01A
+TCGA_BH_A0H5_01A
+TCGA_D8_A27M_01A
+TCGA_AR_A5QP_01A
+TCGA_AO_A12F_01A
+TCGA_EW_A1OV_01A
+TCGA_C8_A1HE_01A
 TCGA_EW_A1PE_01A
-TCGA_EW_A1PF_01A
-TCGA_EW_A1PH_01A
-TCGA_EW_A2FS_01A
-TCGA_EW_A2FW_01A
-TCGA_EW_A6S9_01A
-TCGA_EW_A6SD_01A
-TCGA_GI_A2C9_01A
-TCGA_GM_A2DD_01A
-TCGA_GM_A2DF_01A
-TCGA_GM_A2DH_01A
-TCGA_GM_A2DK_01A
+TCGA_A7_A0DB_01A
+TCGA_E2_A152_01A
+TCGA_B6_A401_01A
+TCGA_D8_A1XW_01A
+TCGA_A8_A07C_01A
+TCGA_EW_A1OW_01A
+TCGA_BH_A0HQ_01A
+TCGA_A2_A04R_01A
+TCGA_E2_A1LG_01A
+TCGA_A2_A0CM_01A
+TCGA_AO_A1KO_01A
+TCGA_B6_A0RH_01A
+TCGA_AO_A03P_01A
+TCGA_BH_A1FU_01A
+TCGA_A1_A0SE_01A
+TCGA_A7_A6VY_01A
+TCGA_AO_A129_01A
+TCGA_E9_A1R0_01A
+TCGA_E2_A158_01A
 TCGA_GM_A2DL_01A
-TCGA_GM_A2DM_01A
-TCGA_GM_A2DN_01A
-TCGA_GM_A3XL_01A
-TCGA_HN_A2NL_01A
-TCGA_HN_A2OB_01A
-TCGA_JL_A3YW_01A
-TCGA_JL_A3YX_01A
-TCGA_LD_A66U_01A
-TCGA_LD_A7W6_01A
-TCGA_LL_A442_01A
-TCGA_LL_A5YN_01A
-TCGA_LL_A5YP_01A
-TCGA_LL_A6FR_01A
-TCGA_MS_A51U_01A
-TCGA_OK_A5Q2_01A
-TCGA_OL_A5D8_01A
-TCGA_OL_A5RW_01A
-TCGA_OL_A5RX_01A
-TCGA_OL_A5RY_01A
-TCGA_OL_A5RZ_01A
-TCGA_OL_A66I_01A
-TCGA_OL_A66O_01A
-TCGA_OL_A6VO_01A
-TCGA_OL_A6VR_01A
-TCGA_PE_A5DE_01A
-TCGA_S3_AA14_01A
-TCGA_S3_AA15_01A
-TCGA_Z7_A8R6_01A
+TCGA_C8_A27B_01A
+TCGA_C8_A27A_01A
+TCGA_E9_A1NE_01A
+TCGA_E9_A1N8_01A
+TCGA_A1_A0SD_01A
+TCGA_E9_A1NI_01A
+TCGA_A2_A0CX_01A
+TCGA_E2_A15D_01A
\ No newline at end of file
diff --git a/tests/simulations_original_BN.npy b/tests/simulations_original_BN.npy
deleted file mode 100644
index bfdefa2..0000000
Binary files a/tests/simulations_original_BN.npy and /dev/null differ
diff --git a/tests/test_breast_cancer.py b/tests/test_breast_cancer.py
index 937878d..1817f62 100644
--- a/tests/test_breast_cancer.py
+++ b/tests/test_breast_cancer.py
@@ -1,5 +1,4 @@
 import os
-import random
 import shutil
 import time
 from typing import List, Optional
@@ -153,13 +152,15 @@ def test_patient_model_simulations(
     # Create patient-specific models
     for patient in TNBC_ID:
         shutil.copytree(path_to_patient("TCGA_3C_AALK_01A"), path_to_patient(f"{patient}"))
+    # Execute patient-specific models
+    simulations = PatientModelSimulations(
+        tests.models.breast.__package__, TNBC_ID if exec_model else [TNBC_ID[0]]
+    )
+    start = time.time()
+    assert simulations.run() is None
+    elapsed = time.time() - start
+    print(f"Computation time for simulating {len(TNBC_ID)} patients: {elapsed/60:.1f} [min]")
     if exec_model:
-        # Execute patient-specific models
-        simulations = PatientModelSimulations(tests.models.breast.__package__, TNBC_ID)
-        start = time.time()
-        assert simulations.run() is None
-        elapsed = time.time() - start
-        print(f"Computation time for simulating {len(TNBC_ID)} patients: {elapsed/60:.1f} [min]")
         # Extract response characteristics and visualize patient classification
         if dynamical_features is None:
             dynamical_features = ["AUC"]
@@ -187,7 +188,7 @@ def test_patient_model_analyses(exec_model: bool = False):
     for patient in TNBC_ID:
         assert os.path.isdir(os.path.join(PATH_TO_MODELS, patient))
     if exec_model:
-        patients = random.sample(TNBC_ID, 1)
+        patients = TNBC_ID
         analyses = PatientModelAnalyses(
             tests.models.breast.__package__,
             patients,
diff --git a/tests/test_error_message.py b/tests/test_error_message.py
deleted file mode 100644
index 18efce3..0000000
--- a/tests/test_error_message.py
+++ /dev/null
@@ -1,57 +0,0 @@
-import os
-import shutil
-
-import pytest
-
-from pasmopy import Text2Model
-from pasmopy.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_fos_model.py b/tests/test_fos_model.py
deleted file mode 100644
index 91480ea..0000000
--- a/tests/test_fos_model.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import os
-import shutil
-
-import numpy as np
-
-from pasmopy 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.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()
-    for obs_name in model.observables:
-        assert os.path.isfile(
-            os.path.join(model.path, "figure", "simulation", "1", f"{obs_name}.png")
-        )
-
-
-def test_cleanup():
-    assert os.path.isdir(PATH_TO_MODEL)
-    shutil.rmtree(PATH_TO_MODEL)
diff --git a/tests/test_model_construction.py b/tests/test_model_construction.py
deleted file mode 100644
index c2dbd6f..0000000
--- a/tests/test_model_construction.py
+++ /dev/null
@@ -1,184 +0,0 @@
-import os
-import shutil
-from typing import Callable
-
-import numpy as np
-import pytest
-
-from pasmopy import Model, Text2Model, run_simulation
-from pasmopy.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.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(n_reaction=2)
-        elif model == "Kholodenko1999":
-            mapk_cascade = Text2Model(
-                os.path.join(
-                    os.path.dirname(__file__),
-                    "text_files",
-                    f"{model}.txt",
-                )
-            )
-            mapk_cascade.to_markdown(n_reaction=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_thermodynamic_restriction.py b/tests/test_thermodynamic_restriction.py
deleted file mode 100644
index 216ecc0..0000000
--- a/tests/test_thermodynamic_restriction.py
+++ /dev/null
@@ -1,54 +0,0 @@
-import os
-import shutil
-
-from pasmopy 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/text_files/Kholodenko1999.txt b/tests/text_files/Kholodenko1999.txt
deleted file mode 100644
index b09eb22..0000000
--- a/tests/text_files/Kholodenko1999.txt
+++ /dev/null
@@ -1,44 +0,0 @@
-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/text_files/abc.txt b/tests/text_files/abc.txt
deleted file mode 100644
index 5b6ca8c..0000000
--- a/tests/text_files/abc.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-A + B <--> AB
-AB + C <--> ABC
-ABC <--> A + BC
-BC <--> B + C
\ No newline at end of file
diff --git a/tests/text_files/duplicate_binding.txt b/tests/text_files/duplicate_binding.txt
deleted file mode 100644
index d5a8397..0000000
--- a/tests/text_files/duplicate_binding.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-A + B <--> AB
-AB <--> A + B
\ No newline at end of file
diff --git a/tests/text_files/fos_model.txt b/tests/text_files/fos_model.txt
deleted file mode 100644
index 493f90e..0000000
--- a/tests/text_files/fos_model.txt
+++ /dev/null
@@ -1,80 +0,0 @@
-@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/text_files/michaelis_menten.txt b/tests/text_files/michaelis_menten.txt
deleted file mode 100644
index af51bc3..0000000
--- a/tests/text_files/michaelis_menten.txt
+++ /dev/null
@@ -1,10 +0,0 @@
-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/text_files/typo_1.txt b/tests/text_files/typo_1.txt
deleted file mode 100644
index 93cc435..0000000
--- a/tests/text_files/typo_1.txt
+++ /dev/null
@@ -1,1 +0,0 @@
-A bindss B <--> AB
\ No newline at end of file
diff --git a/tests/text_files/typo_2.txt b/tests/text_files/typo_2.txt
deleted file mode 100644
index be596aa..0000000
--- a/tests/text_files/typo_2.txt
+++ /dev/null
@@ -1,1 +0,0 @@
-A bnds B <--> AB
\ No newline at end of file