diff --git a/petab/v1/models/__init__.py b/petab/v1/models/__init__.py index a35ad432..bb5dd415 100644 --- a/petab/v1/models/__init__.py +++ b/petab/v1/models/__init__.py @@ -4,12 +4,21 @@ MODEL_TYPE_SBML = "sbml" #: PySB model type as used in a PEtab v2 yaml file as `language`. MODEL_TYPE_PYSB = "pysb" +#: BNGL model type as used in a PEtab v2 yaml file as `language`. +MODEL_TYPE_BNGL = "bngl" known_model_types = { MODEL_TYPE_SBML, MODEL_TYPE_PYSB, + MODEL_TYPE_BNGL, } from .model import Model # noqa F401 -__all__ = ["MODEL_TYPE_SBML", "MODEL_TYPE_PYSB", "known_model_types", "Model"] +__all__ = [ + "MODEL_TYPE_SBML", + "MODEL_TYPE_PYSB", + "MODEL_TYPE_BNGL", + "known_model_types", + "Model", +] diff --git a/petab/v1/models/bngl_model.py b/petab/v1/models/bngl_model.py new file mode 100644 index 00000000..ebdb5003 --- /dev/null +++ b/petab/v1/models/bngl_model.py @@ -0,0 +1,413 @@ +"""Functions for handling BNGL (BioNetGen Language) models. + +``petab`` ships ``sbml`` and ``pysb`` model loaders; this module adds +``bngl`` so that a ``language: bngl`` PEtab problem can be loaded and +validated at the model level (see PEtab-dev/PEtab#436). + +:class:`BnglModel` is a :class:`petab.v1.models.model.Model` backed by a +small, dependency-free BNGL *block reader* (:func:`parse_bngl`). PEtab +validation only ever introspects a model -- it enumerates the model's named +entities -- so the reader does not run BNG2.pl or generate a reaction +network. The single exception is :meth:`BnglModel.is_valid`, which shells +out to ``BNG2.pl --check`` (a parse/semantic check, *without* network +generation) when a BNG2.pl is locatable, and degrades gracefully to +``True`` when no BNG backend is available -- mirroring how the SBML loader +always validates because ``libsbml`` is always present. + +The BNGL entity sets that back the introspection methods were established +against BioNetGen's ``Perl2/`` source (the reference implementation), not +inferred from the PySB analogy: + +* expression symbols are exactly the BNG ``ParamList`` -- parameters, + observables, and global functions; compartments are *not* expression + symbols (they never enter the ``ParamList``); +* the full model-entity namespace additionally includes molecule types, + compartments, and seed species. + +The block scanner is hardened against the BNGL grammar (BioNetGen ``Perl2/``; +cross-checked against ``BNG_vscode_extension`` ``docs/bngl-grammar.md``): line +continuations (a trailing ``\\``), the ``species`` block alias (``begin +species`` = ``begin seed species``), line labels (a numeric index ``1 L0 1`` or +a named ``CD14: ...``), and the seed-species ``$`` clamp marker are honored. It +is kept in sync with PyBNF's sibling reader (``pybnf/petab/_bngl.py``, +ADR-0026) -- the grammar-hardening tests are the anchor against drift. +""" + +from __future__ import annotations + +import os +import re +import shutil +import subprocess +from collections.abc import Iterable +from dataclasses import dataclass +from pathlib import Path + +from ..._utils import _generate_path +from .. import is_valid_identifier +from . import MODEL_TYPE_BNGL +from .model import Model + +__all__ = ["BnglModel", "BnglEntities", "parse_bngl"] + +#: The three keywords that open an observable declaration line. +_OBS_KEYWORDS = frozenset({"Molecules", "Species", "Counter"}) + +#: Short spellings BNG2.pl accepts for a block's canonical (long) name; either +#: spelling opens and closes the same block. The grammar doc +#: (``BNG_vscode_extension/docs/bngl-grammar.md``) also lists ``molecules`` and +#: ``rules``, but BNG2.pl 2.9.3 -- the reference this reader targets -- rejects +#: both ("Could not process block type"), so honoring them would accept models +#: BNG2.pl refuses. Only ``species`` (for ``seed species``, also BNG2.pl's own +#: canonical output spelling) is real. +_BLOCK_ALIASES = { + "seed species": ("species",), +} + + +@dataclass(frozen=True) +class BnglEntities: + """The named entities of a BNGL model that the PEtab layer reads. + + :ivar text: The verbatim BNGL model text. + :ivar parameters: Maps a parameter name to its raw right-hand side -- a + number (``"5"``, ``"6.02e23"``) or an expression (``"2*base"``), + kept verbatim; numeric coercion is the caller's job. + :ivar observable_names: Bare observable names. + :ivar function_names: Bare global-function names (without ``()``). + :ivar molecule_type_names: Bare molecule-type names. + :ivar seed_species: Concrete seed-species pattern strings (verbatim). + :ivar compartment_names: Bare compartment names. + """ + + text: str + parameters: dict[str, str] + observable_names: frozenset[str] + function_names: frozenset[str] + molecule_type_names: frozenset[str] + seed_species: frozenset[str] + compartment_names: frozenset[str] + + +def parse_bngl(text: str) -> BnglEntities: + """Parse BNGL ``text`` into a :class:`BnglEntities`. + + A stdlib ``begin``/``end `` scanner -- no BNG2.pl, no network + generation. Sufficient for PEtab validation, which only introspects the + model's declared entities. + + :param text: The BNGL model text. + :returns: The model's named entities. + """ + parameters: dict[str, str] = {} + for line in _block_lines(text, "parameters"): + name_value = _parameter_name_value(line) + if name_value is not None: + parameters[name_value[0]] = name_value[1] + return BnglEntities( + text=text, + parameters=parameters, + observable_names=_names(text, "observables", _observable_name), + function_names=_names(text, "functions", _function_name), + molecule_type_names=_names( + text, "molecule types", _molecule_type_name + ), + seed_species=_names(text, "seed species", _seed_species_pattern), + compartment_names=_names(text, "compartments", _compartment_name), + ) + + +def _names(text: str, block_name: str, extractor) -> frozenset[str]: + """The non-empty names ``extractor`` yields over a block's lines.""" + return frozenset( + name + for name in ( + extractor(line) for line in _block_lines(text, block_name) + ) + if name + ) + + +def _logical_lines(text: str) -> list[str]: + """The comment-stripped *logical* lines of ``text`` -- physical lines with + BNGL line continuations joined. + + Mirrors BNG2.pl's ``readFile`` (``Perl2/BNGModel.pm``): strip the ``#`` + comment first, then while the line ends with ``\\`` (as its last + non-whitespace character) drop that ``\\`` and append the next + comment-stripped physical line *directly* -- no separating space, so a + token split across the break (``1e\\`` + ``3`` -> ``1e3``) rejoins. + Without this, a continued parameter / function / observable is truncated at + the ``\\`` (e.g. a ``k = \\`` line would read as the value ``"\\"``). + """ + raw_lines = text.splitlines() + out = [] + i, n = 0, len(raw_lines) + while i < n: + line = raw_lines[i].split("#", 1)[0] + i += 1 + while re.search(r"\\\s*$", line): + line = re.sub(r"\\\s*$", "", line) + if i >= n: + break # a dangling continuation at EOF + line += raw_lines[i].split("#", 1)[0] + i += 1 + out.append(line.strip()) + return out + + +def _block_lines(text: str, block_name: str) -> list[str]: + """The comment-stripped, non-blank lines inside ``begin``/``end``. + + ``block_name`` is the canonical (long) spelling; a BNG2.pl-accepted alias + for it (only ``species`` for ``seed species``; see :data:`_BLOCK_ALIASES`) + opens and closes the same block. Lines are logical lines -- continuations + joined (see :func:`_logical_lines`). + """ + names = "|".join( + re.escape(name) + for name in (block_name, *_BLOCK_ALIASES.get(block_name, ())) + ) + begin = re.compile(rf"^begin\s+(?:{names})\b", re.IGNORECASE) + end = re.compile(rf"^end\s+(?:{names})\b", re.IGNORECASE) + lines = [] + in_block = False + for line in _logical_lines(text): + if begin.match(line): + in_block = True + elif end.match(line): + in_block = False + elif in_block and line: + lines.append(line) + return lines + + +def _strip_line_label(line: str) -> str: + """Drop a leading BNGL line label so the entity, not the label, is read. + + ``LineLabel = {Digit}, WS | Name, ":", [WS]`` (grammar) -- either a numeric + index (the legacy ``.net``-style ``1 L0 1`` form) or a named label + (``CD14: CD14(...)``). A valid BNGL identifier starts with a letter, so a + leading digit-run is always an index; a compartment prefix is ``@Name:`` + (with the ``@``), so a bare ``Name:`` at line start is unambiguously a + label. + """ + match = re.match(r"^\d+\s+(.*)$", line) or re.match( + r"^[A-Za-z]\w*:\s+(.*)$", line + ) + return match.group(1) if match else line + + +def _parameter_name_value(line: str) -> tuple[str, str] | None: + """``(name, rhs)`` for a ``[LineLabel] Name (WS|"=") MathExpr`` line.""" + line = _strip_line_label(line) + match = re.match(r"^(\w+)\s*=\s*(.+)$", line) or re.match( + r"^(\w+)\s+(.+)$", line + ) + return (match.group(1), match.group(2).strip()) if match else None + + +def _observable_name(line: str) -> str | None: + """The name in a `` `` observable line.""" + tokens = line.split() + if len(tokens) >= 2 and tokens[0] in _OBS_KEYWORDS: + return tokens[1] + return None + + +def _function_name(line: str) -> str | None: + """The name in a ``() = ...`` global-function line.""" + match = re.match(r"(\w+)\s*\(", line) or re.match(r"(\w+)\s*=", line) + return match.group(1) if match else None + + +def _molecule_type_name(line: str) -> str | None: + """The name in a ``(...)`` molecule-type line.""" + match = re.match(r"(\w+)", line) + return match.group(1) if match else None + + +def _seed_species_pattern(line: str) -> str | None: + """The species pattern in a ``[LineLabel] ["$"] `` line. + + A leading line label (numeric index ``1 A() 100`` or named + ``CD14: CD14(...)``; see :func:`_strip_line_label`) is dropped first so the + label is not mistaken for the species. A leading ``$`` (the fixed/clamped- + concentration marker, ``SeedSpeciesDefn = ["$"], Species, WS, + MathExpression``) is a modifier, not part of the species identity, so it + too is stripped: ``$counter() 10`` + enumerates the state variable ``counter()``. + """ + line = _strip_line_label(line) + if line.startswith("$"): + line = line[1:].lstrip() + tokens = line.split() + return tokens[0] if tokens else None + + +def _compartment_name(line: str) -> str | None: + """The name in a `` [outside]`` line.""" + tokens = line.split() + return tokens[0] if tokens else None + + +class BnglModel(Model): + """PEtab wrapper for BNGL models (introspection only; no simulation).""" + + type_id = MODEL_TYPE_BNGL + + def __init__( + self, + model: BnglEntities, + model_id: str | None = None, + rel_path: Path | str | None = None, + base_path: str | Path | None = None, + ): + super().__init__() + + self.rel_path = rel_path + self.base_path = base_path + + self.model = model + self._model_id = model_id + + if not is_valid_identifier(self._model_id): + raise ValueError( + f"Model ID '{self._model_id}' is not a valid identifier. " + "Either provide a valid identifier or rename the model file " + "to a valid PEtab model identifier." + ) + + @staticmethod + def from_file( + filepath_or_buffer, + model_id: str | None = None, + base_path: str | Path | None = None, + ) -> BnglModel: + path = Path(_generate_path(filepath_or_buffer, base_path)) + text = path.read_text(encoding="utf-8", errors="replace") + return BnglModel( + model=parse_bngl(text), + model_id=model_id or path.stem, + rel_path=filepath_or_buffer, + base_path=base_path, + ) + + def to_file(self, filename: str | Path | None = None) -> None: + target = filename or _generate_path(self.rel_path, self.base_path) + with open(target, "w", encoding="utf-8") as f: + f.write(self.model.text) + + @property + def model_id(self): + return self._model_id + + @model_id.setter + def model_id(self, model_id): + self._model_id = model_id + + # -- parameters --------------------------------------------------------- + + def get_parameter_ids(self) -> Iterable[str]: + return list(self.model.parameters) + + def get_parameter_value(self, id_: str) -> float: + try: + rhs = self.model.parameters[id_] + except KeyError as e: + raise ValueError(f"Parameter {id_} does not exist.") from e + try: + return float(rhs) + except ValueError as e: + raise NotImplementedError( + f"Parameter '{id_}' has an expression value '{rhs}'. " + "Evaluating a BNGL parameter expression requires BNG2.pl / " + "network generation, which is out of scope for the " + "introspection-only BnglModel." + ) from e + + def get_free_parameter_ids_with_values( + self, + ) -> Iterable[tuple[str, float]]: + out = [] + for name, rhs in self.model.parameters.items(): + try: + out.append((name, float(rhs))) + except ValueError: + # An expression-valued parameter has no introspection-grade + # value; skip it rather than evaluate the expression. + continue + return out + + def get_valid_parameters_for_parameter_table(self) -> Iterable[str]: + # All parameters are allowed in the parameter table. + return list(self.model.parameters) + + # -- model-entity namespaces (verified against BioNetGen) --------------- + + def has_entity_with_id(self, entity_id) -> bool: + # The full declared-identifier namespace. + return ( + entity_id in self.model.parameters + or entity_id in self.model.observable_names + or entity_id in self.model.function_names + or entity_id in self.model.molecule_type_names + or entity_id in self.model.compartment_names + or entity_id in self.model.seed_species + ) + + def get_valid_ids_for_condition_table(self) -> Iterable[str]: + return list(self.model.parameters) + list(self.model.compartment_names) + + def symbol_allowed_in_observable_formula(self, id_: str) -> bool: + # The BNG ParamList: parameters, observables, global functions only. + return ( + id_ in self.model.parameters + or id_ in self.model.observable_names + or id_ in self.model.function_names + ) + + def is_state_variable(self, id_: str) -> bool: + # At introspection grade only the concrete seed species are known; + # the full species set is a network-generation product. + return id_ in self.model.seed_species + + # -- validity ----------------------------------------------------------- + + def is_valid(self) -> bool: + # Real BNG2.pl --check (parse/semantic validation, no network + # generation) when locatable, else True -- never a false failure + # where no BNG backend is available. + bng2 = _locate_bng2() + if bng2 is None or self.rel_path is None: + return True + path = Path(_generate_path(self.rel_path, self.base_path)) + if not path.is_file(): + # No local model file to check (e.g. a buffer-loaded model). + return True + # Resolve to an absolute path so BNG2.pl finds the model regardless + # of the working directory it runs in (its output stays next to the + # model via ``cwd``). + path = path.resolve() + try: + result = subprocess.run( # noqa: S603 + [bng2, "--check", str(path)], + capture_output=True, + text=True, + timeout=120, + cwd=str(path.parent), + ) + except (OSError, subprocess.SubprocessError): + # A tooling hiccup must not masquerade as an invalid model. + return True + return result.returncode == 0 + + +def _locate_bng2() -> str | None: + """A path to ``BNG2.pl`` via ``BNGPATH`` or ``PATH``, else ``None``.""" + bngpath = os.environ.get("BNGPATH") + if bngpath: + candidate = Path(bngpath) / "BNG2.pl" + if candidate.is_file(): + return str(candidate) + return shutil.which("BNG2.pl") diff --git a/petab/v1/models/model.py b/petab/v1/models/model.py index 96613757..d1a9da9d 100644 --- a/petab/v1/models/model.py +++ b/petab/v1/models/model.py @@ -150,7 +150,12 @@ def model_factory( :param model_id: PEtab model ID for the given model :returns: A :py:class:`Model` instance representing the given model """ - from . import MODEL_TYPE_PYSB, MODEL_TYPE_SBML, known_model_types + from . import ( + MODEL_TYPE_BNGL, + MODEL_TYPE_PYSB, + MODEL_TYPE_SBML, + known_model_types, + ) if model_language == MODEL_TYPE_SBML: from .sbml_model import SbmlModel @@ -166,6 +171,13 @@ def model_factory( filepath_or_buffer, model_id=model_id, base_path=base_path ) + if model_language == MODEL_TYPE_BNGL: + from .bngl_model import BnglModel + + return BnglModel.from_file( + filepath_or_buffer, model_id=model_id, base_path=base_path + ) + if model_language in known_model_types: raise NotImplementedError( f"Unsupported model format: {model_language}" diff --git a/petab/v2/models/bngl_model.py b/petab/v2/models/bngl_model.py new file mode 100644 index 00000000..d4a291ab --- /dev/null +++ b/petab/v2/models/bngl_model.py @@ -0,0 +1,3 @@ +"""Functions for handling BNGL models""" + +from ...v1.models.bngl_model import * # noqa: F401, F403 diff --git a/tests/v1/bngl/parabola.bngl b/tests/v1/bngl/parabola.bngl new file mode 100644 index 00000000..aa1e0245 --- /dev/null +++ b/tests/v1/bngl/parabola.bngl @@ -0,0 +1,30 @@ +begin model + + begin parameters + v1 5 + v2 5 + v3 5 + k_cond 1 + end parameters + + begin molecule types + counter() + end molecule types + + begin seed species + counter() -10 + end seed species + + begin observables + Molecules x counter() + end observables + + begin functions + y()=v1*(x^2)+(v2*x)+v3 + end functions + + begin reaction rules + 0->counter() 1 + end reaction rules + +end model diff --git a/tests/v1/bngl_corpus/An_2009.bngl b/tests/v1/bngl_corpus/An_2009.bngl new file mode 100644 index 00000000..fd28e17c --- /dev/null +++ b/tests/v1/bngl_corpus/An_2009.bngl @@ -0,0 +1,248 @@ +begin parameters +LPS_MD2_Bind 0.001 +LPS_MD2_Unbind 0.1 +LPS_CD14_Bind 0.001 +LPS_CD14_Unbind 0.1 +CD14_MD2_Bind 0.001 +CD14_MD2_Unbind 0.1 +LPS_TLR4_Bind 0.001 +LPS_TLR4_Unbind 0.1 +CD14_TLR4_Bind 0.001 +CD14_TLR4_Unbind 0.1 +MD2_TLR4_Bind 0.001 +MD2_TLR4_Unbind 0.1 +TLR4_Complex_Dimer_Bind 0.001 +TLR4_Complex_Dimer_Unbind 0.1 +RP1_TRIF_Bind 0.001 +RP1_TRIF_Unbind 0.1 +TRIF_TRAF6_Bind 0.001 +TRIF_TRAF6_Unbind 0.1 +RP1_TRAF6_Bind 0.001 +RP1_TRAF6_Unbind 0.1 +TLR4_TRAM_Bind 0.001 +TLR4_TRAM_Unbind 0.1 +TLR4TRAM_TRIF_Bind 0.001 +TLR4TRAM_TRIF_Unbind 0.1 +TRAF6_TRIF_Bind 0.001 +TRAF6_TRIF_Unbind 0.1 +TRAF6TRIF_TAK1_Activate 0.001 +MyD88_IRAK4_Bind 0.001 +MyD88_IRAK4_Unbind 0.1 +MyD88_IRAK1_Bind 0.001 +MyD88_IRAK1_Unbind 0.1 +IRAK1_IRAK4_Bind 0.001 +IRAK1_IRAK4_Unbind 0.1 +TLR4_MAL_Bind 0.001 +TLR4_MAL_Unbind 0.1 +TLR4MAL_MyD88_Bind 0.001 +TLR4MAL_MyD88_Unbind 0.1 +MyD88IRAK1_TRAF6_Bind 0.001 +MyD88IRAK1_TRAF6_B_Unbind 0.1 +MyD88IRAK1TRAF6_TAK1_Activate 0.001 +TRAF6_MyD88IRAK1_Bind 0.001 +TRAF6_MyD88IRAK1_Unbind 0.1 +TAK1_Ikk_Complex_Activate 0.001 +Ikk_complex_IkB_Phos 0.00001 +IkB_Proteasome23_Degrade 0.1 +p65_p50_Bind 0.001 +p65_p50_Unbind 0.1 +NFkB_DNA_A20_Bind 0.001 +A20_Transcription_Execute 1 +A20_TRAF6_Bind 0.001 +A20_TRAF6_Unbind 0.1 +NFkB_DNA_A20_Unbind 0.1 +NFkB_DNA_TNF_Bind 0.001 +NFkB_DNA_TNF_Unbind 0.1 +TNF_Transcription_Execute 1 +CD14_Init 10000 +MD2_Init 10000 +TLR_Init 10000 +LPS_Init 100 +TRAM_Init 10000 +MAL_Init 10000 +TRIF_Init 10000 +MyD88_Init 10000 +RP1_Init 10000 +IRAK1_Init 10000 +IRAK4_Init 10000 +TRAF6_Init 10000 +TAK1_Init 10000 +Ikk_Complex_Init 10000 +Proteasome23_Init 10000 +p65_Init 10000 +IkB_Init 10000 +p50_Init 10000 +DNA 2 +A20_Translation_Execute 0.1 +TNF_Translation_Execute 0.1 +TAK1_Degradation 0 +NFkB_Degredation 0 +A20_MyD88IRAK1TRAF6_Degrade 10 +A20_TRAF6TRIFRP1_Degrade 10 +A20_Init 0 +A20_Preconditioned 0 +Ikk_Degradation_Rate 0 +NFkB_DNA_IkB_Bind 0.001 +NFkB_DNA_IkB_Unbind 0.01 +IkB_Transcription_Execute 1 +IkB_Translation_Execute 0.1 +A20_IkkAct_Deactivate 10 +IkB_DegradeNFkB 0.001 +NFkB_Inactive_Cytoplasm 10000 +NFkB_IkB_Bind 0.001 +NFkB_Translocation_Nucleus 0.01 +NFkB_IkB_Unbind 0.0001 +TAK1_Deactivation 0.10 +Ikk_Deactivation 0 +TNF_Degrade 0.0005 +A20_Degrade 0.0001 +end parameters + +begin molecule types + CD14(TLR4,MD2,LPS) + MD2(CD14,TLR4,LPS) + TLR4(MAL,TRAM,TLR4,CD14,MD2,LPS) + TRAM(TLR4,TRIF) + TRIF(TRAM,TRAF6,RIP1,TRAF4,SARM) + SARM(TRIF) + TRAF4(TRAF6,TAK1,TRIF) + IRAK1(IRAK4,MyD88,Tollip,TRAF6) + Tollip(IRAK1) + IRAK4(Myd88,IRAKM,IRAK1) + IRAKM(IRAK4) + RP1(TRIF,TRAF6,TAK1,p38) + TRAF6(IRAK1,TRIF,RP1,TAK1,TRAF4,A20,JNK,p38) + A20(TRAF6) + MyD88(MAL,IRAK1,IRAK4,MyD88s) + MyD88s(MyD88,IRAK1) + MAL(TLR4,MyD88,SOCS1) + LPS(MD2,TLR4,CD14,LPS) + TAK1(TRAF6,Activation~No~Yes) + Ikk_Complex(Activation~No~Yes) + IkB(Phos~No~Yes,p65,p50,Degrade~No~Yes) + Proteasome26s(IkB) + TNF(TNFr) + TNFmRNA(Translation~On~Off) + A20mRNA(Translation~On~Off) + iNOSmRNA(Translation~On~Off) + IkBmRNA(Translation~On~Off) + DNA(A20,TNF,iNOS,IL10,IkB,c,c) + Trash(c) + Administer(c) + NFkB(Transcription~No~Yes,Activation~No~Yes,Location~Cytoplasm~Nucleus) +end molecule types + +begin seed species +CD14: CD14(TLR4,MD2,LPS) CD14_Init +MD2: MD2(CD14,TLR4,LPS) MD2_Init +TLR4: TLR4(MAL,TRAM,TLR4,CD14,MD2,LPS) TLR_Init +TRAM: TRAM(TLR4,TRIF) TRAM_Init +MAL: MAL(TLR4,MyD88,SOCS1) MAL_Init +TRIF: TRIF(TRAM,TRAF6,RIP1,SARM,TRAF4) TRIF_Init +MyD88: MyD88(MAL,IRAK1,IRAK4,MyD88s) MyD88_Init +RP1: RP1(TRIF,TRAF6,TAK1,p38) RP1_Init +IRAK1: IRAK1(IRAK4,MyD88,Tollip,TRAF6) IRAK1_Init +IRAK4: IRAK4(Myd88,IRAKM,IRAK1) IRAK4_Init +TRAF6: TRAF6(IRAK1,TRIF,RP1,TAK1,TRAF4,A20,JNK,p38) TRAF6_Init +TAK1: TAK1(TRAF6,Activation~No) TAK1_Init +Ikk_Complex: Ikk_Complex(Activation~No) Ikk_Complex_Init +Proteasome23: Proteasome26s(IkB) Proteasome23_Init +IkB: IkB(Phos~No,p65,p50,Degrade~No) IkB_Init +DNA: DNA(A20,TNF,iNOS,IL10,IkB,c,c) DNA +LPS: LPS(MD2,TLR4,CD14,LPS) LPS_Init +A20: A20(TRAF6) A20_Preconditioned +NFkB_Inactive_Cytoplasm: NFkB(Transcription~No,Activation~No,Location~Cytoplasm) NFkB_Inactive_Cytoplasm +end seed species + +begin reaction rules +LPS_MD2: LPS(MD2,TLR4,CD14,LPS)+MD2(CD14,TLR4,LPS)<->LPS(MD2!0,TLR4,CD14,LPS).MD2(CD14,TLR4,LPS!0) LPS_MD2_Bind,LPS_MD2_Unbind +LPS_MD2_CD14: LPS(MD2!0,TLR4,CD14,LPS).MD2(CD14,TLR4,LPS!0)+CD14(TLR4,MD2,LPS)<->LPS(MD2!0,TLR4,CD14!1,LPS).MD2(CD14,TLR4,LPS!0).CD14(TLR4,MD2,LPS!1) LPS_CD14_Bind,LPS_CD14_Unbind +LPS_MD2_CD14_TLR4: LPS(MD2!0,TLR4,CD14!1,LPS).MD2(CD14,TLR4,LPS!0).CD14(TLR4,MD2,LPS!1) + TLR4(MAL,TRAM,TLR4,CD14,MD2,LPS) <-> LPS(MD2!0,TLR4!2,CD14!1,LPS).MD2(CD14!3,TLR4!4,LPS!0).CD14(TLR4!5,MD2!3,LPS!1).TLR4(MAL,TRAM,TLR4,CD14!5,MD2!4,LPS!2) LPS_TLR4_Bind,LPS_TLR4_Unbind +TLR4_Dimerization: TLR4(TLR4,CD14!+,LPS!+,MD2!+)+TLR4(TLR4,CD14!+,MD2!+,LPS!+)<->TLR4(TLR4!0,CD14!+,LPS!+,MD2!+).TLR4(TLR4!0,CD14!+,LPS!+,MD2!+) TLR4_Complex_Dimer_Bind,TLR4_Complex_Dimer_Unbind +TLR4_TRAM: TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL)+TRAM(TLR4,TRIF)<->TRAM(TLR4!0,TRIF).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM!0,MAL) TLR4_TRAM_Bind,TLR4_TRAM_Unbind +TLR4TRAM_TRIF: TRAM(TLR4!0,TRIF).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM!0,MAL)+TRIF(TRAM,TRAF6!+,RIP1!+,TRAF4,SARM)<->TRAM(TLR4!0,TRIF!1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM!0,MAL).TRIF(TRAM!1,TRAF6!+,RIP1!+,TRAF4,SARM) TLR4TRAM_TRIF_Bind,TLR4TRAM_TRIF_Unbind +TLR4TRAMTRIFTRAF6_TAK1: TRAM(TLR4!0,TRIF!1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM!0,MAL).TRIF(TRAM!1,TRAF6!+,RIP1!+,TRAF4,SARM)+TAK1(TRAF6,Activation~No)->TRAM(TLR4!0,TRIF!1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM!0,MAL).TRIF(TRAM!1,TRAF6!+,RIP1!+,TRAF4,SARM)+TAK1(TRAF6,Activation~Yes) TRAF6TRIF_TAK1_Activate +MyD88_IRAK4: MyD88(MAL,IRAK1,IRAK4,MyD88s)+IRAK4(Myd88,IRAKM,IRAK1)<->MyD88(MAL,IRAK1,IRAK4!0,MyD88s).IRAK4(Myd88!0,IRAKM,IRAK1) MyD88_IRAK4_Bind,MyD88_IRAK4_Unbind +MyD88_IRAK1: MyD88(MAL,IRAK1,IRAK4!0,MyD88s).IRAK4(Myd88!0,IRAKM,IRAK1)+IRAK1(IRAK4,MyD88,Tollip,TRAF6)<->IRAK1(IRAK4,MyD88!0,Tollip,TRAF6).MyD88(MAL,IRAK1!0,IRAK4!1,MyD88s).IRAK4(Myd88!1,IRAKM,IRAK1) MyD88_IRAK1_Bind,MyD88_IRAK1_Unbind +TNF_Translation: TNFmRNA(Translation~On)->TNFmRNA(Translation~Off)+TNF(TNFr) TNF_Translation_Execute +A20_Translation: A20mRNA(Translation~On)->A20mRNA(Translation~Off)+A20(TRAF6) A20_Translation_Execute +TAK1_Degrade: TAK1(TRAF6,Activation~Yes)->Trash(c) TAK1_Degradation +TLR4_MAL: TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL)+MAL(TLR4,MyD88,SOCS1)<->TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!0).MAL(SOCS1,TLR4!0,MyD88) TLR4_MAL_Bind,TLR4_MAL_Unbind + +TLR4MAL_MyD88: TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!0).MAL(TLR4!0,MyD88,SOCS1) + IRAK1(IRAK4,MyD88!2,Tollip,TRAF6).MyD88(MAL,IRAK1!2,IRAK4!1,MyD88s).IRAK4(Myd88!1,IRAKM,IRAK1) <-> TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!0).MAL(TLR4!0,MyD88!4,SOCS1).IRAK1(IRAK4,MyD88!2,Tollip,TRAF6).MyD88(MAL!4,IRAK1!2,IRAK4!3,MyD88s).IRAK4(Myd88!3,IRAKM,IRAK1) TLR4MAL_MyD88_Bind,TLR4MAL_MyD88_Unbind +MyD88IRAK1_TRAF6: TRAF6(IRAK1,TRIF,RP1,TRAF4,A20,JNK,p38,TAK1)+TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!0).MAL(TLR4!0,MyD88!1,SOCS1).IRAK1(IRAK4,MyD88!2,Tollip,TRAF6).MyD88(MAL!1,IRAK1!2,IRAK4!3,MyD88s).IRAK4(Myd88!3,IRAKM,IRAK1)<->TRAF6(IRAK1!4,TRIF,RP1,TRAF4,A20,JNK,p38,TAK1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!0).MAL(TLR4!0,MyD88!1,SOCS1).IRAK1(IRAK4,MyD88!2,Tollip,TRAF6!4).MyD88(MAL!1,IRAK1!2,IRAK4!3,MyD88s).IRAK4(Myd88!3,IRAKM,IRAK1) MyD88IRAK1_TRAF6_Bind,MyD88IRAK1_TRAF6_B_Unbind + +A20_MyD88IRAK1TRAF6: TRAF6(IRAK1!0,TRIF,RP1,TRAF4,A20,JNK,p38,TAK1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!1).MAL(TLR4!1,MyD88!2,SOCS1).IRAK1(IRAK4,MyD88!3,Tollip,TRAF6!0).MyD88(MAL!2,IRAK1!3,IRAK4!4,MyD88s).IRAK4(Myd88!4,IRAKM,IRAK1)+A20(TRAF6)->TRAF6(IRAK1,TRIF,RP1,TAK1,TRAF4,A20,JNK,p38)+TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!1).MAL(TLR4!1,MyD88!2,SOCS1).IRAK1(IRAK4,MyD88!3,Tollip,TRAF6).MyD88(MAL!2,IRAK1!3,IRAK4!4,MyD88s).IRAK4(Myd88!4,IRAKM,IRAK1)+A20(TRAF6) A20_MyD88IRAK1TRAF6_Degrade +MyD88IRAK1TRAF6_TAK1: TRAF6(IRAK1!0,TRIF,RP1,TRAF4,A20,JNK,p38,TAK1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!1).MAL(TLR4!1,MyD88!2,SOCS1).IRAK1(IRAK4,MyD88!3,Tollip,TRAF6!0).MyD88(MAL!2,IRAK1!3,IRAK4!4,MyD88s).IRAK4(Myd88!4,IRAKM,IRAK1)+TAK1(TRAF6,Activation~No)->TRAF6(IRAK1!0,TRIF,RP1,TRAF4,A20,JNK,p38,TAK1).TLR4(TLR4!+,CD14!+,LPS!+,MD2!+,TRAM,MAL!1).MAL(TLR4!1,MyD88!2,SOCS1).IRAK1(IRAK4,MyD88!3,Tollip,TRAF6!0).MyD88(MAL!2,IRAK1!3,IRAK4!4,MyD88s).IRAK4(Myd88!4,IRAKM,IRAK1)+TAK1(TRAF6,Activation~Yes) MyD88IRAK1TRAF6_TAK1_Activate +IkB_Translation: IkBmRNA(Translation~On)->IkBmRNA(Translation~Off)+IkB(Phos~No,p65,p50,Degrade~No) IkB_Translation_Execute +TRIF_RP1: TRIF(TRAM,TRAF6,RIP1,TRAF4,SARM)+RP1(TRIF,TRAF6,TAK1,p38)<->RP1(TRIF!0,TRAF6,TAK1,p38).TRIF(TRAM,TRAF6,RIP1!0,TRAF4,SARM) RP1_TRIF_Bind,RP1_TRIF_Unbind +TRIF_TRAF6: RP1(TRIF!0,TRAF6,TAK1,p38).TRIF(TRAM,TRAF6,RIP1!0,TRAF4,SARM)+TRAF6(IRAK1,TRIF,RP1,TAK1,TRAF4,A20,JNK,p38)<->RP1(TRIF!0,TRAF6!1,TAK1,p38).TRIF(TRAM,TRAF6!2,RIP1!0,TRAF4,SARM).TRAF6(IRAK1,TRIF!2,RP1!1,TAK1,TRAF4,A20,JNK,p38) TRIF_TRAF6_Bind,TRIF_TRAF6_Unbind +A20_IkkAct_Deactivate: A20(TRAF6)+Ikk_Complex(Activation~Yes)->A20(TRAF6)+Ikk_Complex(Activation~No) A20_IkkAct_Deactivate +A20_TRAF6TRIFRP1_Degrade: A20(TRAF6)+RP1(TRIF!0,TRAF6!1,TAK1,p38).TRIF(TRAM,TRAF6!2,RIP1!0,TRAF4,SARM).TRAF6(IRAK1,TRIF!2,RP1!1,TAK1,TRAF4,A20,JNK,p38)->TRAF6(IRAK1,TRIF,RP1,TAK1,TRAF4,A20,JNK,p38)+RP1(TRIF,TRAF6,TAK1,p38)+TRIF(TRAM,TRAF6,RIP1,TRAF4,SARM)+A20(TRAF6) A20_TRAF6TRIFRP1_Degrade + +Ikk_complex_IkB_Phos: Ikk_Complex(Activation~Yes)+NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No)->Ikk_Complex(Activation~Yes)+NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm).IkB(p65!0,p50!1,Degrade~No,Phos~Yes) Ikk_complex_IkB_Phos + +NFkB_Translocation_Nucleus: NFkB(Transcription~No,Activation~Yes,Location~Cytoplasm)<->NFkB(Transcription~No,Activation~Yes,Location~Nucleus) NFkB_Translocation_Nucleus,NFkB_Translocation_Nucleus +NFkB_DNA_A20: NFkB(Transcription~No,Activation~Yes,Location~Nucleus)+DNA(A20)<->DNA(A20!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) NFkB_DNA_A20_Bind,NFkB_DNA_A20_Unbind + +IkB_DegradeNFkBA20: DNA(A20!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)+IkB(Phos~No,p65,p50,Degrade~No)->DNA(A20)+NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No) IkB_DegradeNFkB + +NFkB_DNA_TNF: NFkB(Transcription~No,Activation~Yes,Location~Nucleus)+DNA(TNF)<->DNA(TNF!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) NFkB_DNA_TNF_Bind,NFkB_DNA_TNF_Unbind +NFkB_DNA_IkB: NFkB(Transcription~No,Activation~Yes,Location~Nucleus)+DNA(IkB)<->DNA(IkB!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) NFkB_DNA_IkB_Bind,NFkB_DNA_IkB_Unbind +IkB_Proteasome23_Degrade: Proteasome26s(IkB)+NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm).IkB(Phos~Yes,p65!0,p50!1,Degrade~No)->IkB(Phos~Yes,p65,p50,Degrade~Yes!0).Proteasome26s(IkB!0)+NFkB(Transcription~No,Activation~Yes,Location~Cytoplasm) IkB_Proteasome23_Degrade +#NFkB_IkB_Bind: NFkB(Location~Cytoplasm,Activation~*)+IkB(Phos~No,p65,p50,Degrade~No)<->NFkB(Activation~*!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No) NFkB_IkB_Bind,NFkB_IkB_Unbind + +NFkB_IkB_Bind: NFkB(Location~Cytoplasm,Activation)+IkB(Phos~No,p65,p50,Degrade~No)<->NFkB(Activation!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No) NFkB_IkB_Bind,NFkB_IkB_Unbind + +Proteasome23_Release: IkB(Degrade~Yes!0).Proteasome26s(IkB!0)->IkB(Degrade~Yes)+Proteasome26s(IkB) IkB_Proteasome23_Degrade +IkB_Transcription_Execute: DNA(IkB!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)->IkBmRNA(Translation~On)+DNA(IkB!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) IkB_Transcription_Execute +A20_Transcription_Execute: DNA(A20!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)->A20mRNA(Translation~On)+DNA(A20!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) A20_Transcription_Execute +TNF_Transcription_Execute: DNA(TNF!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)->TNFmRNA(Translation~On)+DNA(TNF!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) TNF_Transcription_Execute +TAK1_Deactivation: TAK1(TRAF6,Activation~Yes)->TAK1(TRAF6,Activation~No) TAK1_Deactivation +Ikk_Deactivation: Ikk_Complex(Activation~Yes)->Ikk_Complex(Activation~No) Ikk_Deactivation +TAK1_Ikk_Complex_Activate: TAK1(TRAF6,Activation~Yes)+Ikk_Complex(Activation~No)->TAK1(TRAF6,Activation~Yes)+Ikk_Complex(Activation~Yes) TAK1_Ikk_Complex_Activate +TNF_Degrade: TNF(TNFr)->Trash(c) TNF_Degrade +A20_Degrade: A20(TRAF6)->Trash(c) A20_Degrade +IkB_DegradeNFkBDNAIkB: DNA(IkB!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)+IkB(Phos~No,p65,p50,Degrade~No)->DNA(IkB)+IkB(Phos~No,p65!0,p50!1,Degrade~No).NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm) IkB_DegradeNFkB +IkB_DegradeNFkBDNA_TNF: DNA(TNF!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus)+IkB(Phos~No,p65,p50,Degrade~No)->DNA(TNF)+IkB(Phos~No,p65!0,p50!1,Degrade~No).NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm) IkB_DegradeNFkB +end reaction rules + +begin observables +Molecules TNF TNF(TNFr) +Molecules Activated_TAK1 TAK1(TRAF6,Activation~Yes) +Molecules Activated_Ikk_complex Ikk_Complex(Activation~Yes) +Molecules A20 A20(TRAF6) +Molecules NFkB_Active_Cyto NFkB(Transcription~No,Activation~Yes,Location~Cytoplasm) +Molecules NFkB_Active_Nucleus NFkB(Transcription~No,Activation~Yes,Location~Nucleus) +Molecules IkB_Degraded IkB(Degrade~Yes) +Molecules IkB_active IkB(Degrade~No) +#Molecules NFkB_Inactive NFkB(Activation~*!+) + +Molecules NFkB_Inactive NFkB(Activation!+) +Molecules NonBoundNonPhos_IkB IkB(Phos~No,p65,p50,Degrade~No) +Molecules IkBmRNA_Off IkBmRNA(Translation~Off) +Molecules NFkB_DNA_IkB DNA(IkB!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) +Molecules Phos_IkB_NFkB NFkB(Transcription~No,Activation~No!0!1,Location~Cytoplasm).IkB(Phos~Yes,p65!0,p50!1,Degrade~No) +Molecules IkB_Prot26s IkB(Phos~Yes,p65,p50,Degrade~Yes!0).Proteasome26s(IkB!0) + +#Molecules Unbound_Cyto_NFkB NFkB(Location~Cytoplasm,Activation) +Molecules Unbound_Cyto_NFkB NFkB(Location~Cytoplasm,Activation) + +#Molecules Inactive_Cyto_NFkB NFkB(Activation~*!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No) +Molecules Inactive_Cyto_NFkB NFkB(Activation!0!1,Location~Cytoplasm).IkB(Phos~No,p65!0,p50!1,Degrade~No) + +Molecules TNFmRNA_Off TNFmRNA(Translation~Off) +Molecules TNF_NFkB_DNA DNA(TNF!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) +Molecules A20_NFkB_DNA DNA(A20!0).NFkB(Transcription~Yes!0,Activation~Yes,Location~Nucleus) +end observables + +generate_network({overwrite=>1}); + +setConcentration("LPS(MD2,TLR4,CD14,LPS)",0); +simulate_ode({suffix=>"equil",t_end=>50000,n_steps=>10,atol=>1e-12,rtol=>1e-12,sparse=>1,steady_state=>1}); + +#BEGIN SIMULATION +setConcentration("LPS(MD2,TLR4,CD14,LPS)","LPS_Init"); +writeSBML(); +writeMfile(); +simulate_ode({t_end=>100000,n_steps=>500,atol=>1e-12,rtol=>1e-12,sparse=>0}); diff --git a/tests/v1/bngl_corpus/Barua_2009__PATCHED.bngl b/tests/v1/bngl_corpus/Barua_2009__PATCHED.bngl new file mode 100644 index 00000000..8b49d708 --- /dev/null +++ b/tests/v1/bngl_corpus/Barua_2009__PATCHED.bngl @@ -0,0 +1,52 @@ +begin parameters +kon_dimer 1.0 # units = 1/(uM.s) +koff_dimer 0.1 # units = 1/s +kon_SH2 1.0 # units = 1/(uM.s) +koff_SH2 0.1 # units = 1/s +kphos_slow 0.1 # units = 1/s +kphos_fast 1.0 # units = 1/s +Jtot 0.000014 # units = uM +Stot 0.1 # units = uM +end parameters + +begin molecule types +S(SH2,DD) + +J(Y1~P,Y~U~P) + +end molecule types + +begin species +S(SH2,DD) Stot +J(Y1~P,Y~U) Jtot + +end species + +begin reaction rules + +# JAK2-SH2B interaction +J(Y1~P) + S(SH2) <-> J(Y1~P!1).S(SH2!1) kon_SH2, koff_SH2 + +# SH2B dimerization +S(DD) + S(DD) <-> S(DD!1).S(DD!1) kon_dimer, koff_dimer + +# JAK2 phosphorylation +J(Y~U,Y1!1).S(SH2!1,DD!2).S(DD!2,SH2!3).J(Y1!3,Y~U) -> J(Y~P,Y1!1).S(SH2!1,DD!2).S(DD!2,SH2!3).J(Y1!3,Y~U) kphos_slow + +J(Y~U,Y1!1).S(SH2!1,DD!2).S(DD!2,SH2!3).J(Y1!3,Y~P) -> J(Y~P,Y1!1).S(SH2!1,DD!2).S(DD!2,SH2!3).J(Y1!3,Y~P) kphos_fast + + +end reaction rules + +begin observables + Molecules J_mono J(Y1~P) + Molecules JS J(Y1~P!1).S(SH2!1,DD) + Molecules JSS J(Y1~P!1).S(SH2!1,DD!2).S(SH2,DD!2) + Molecules JSSJ J.J + Molecules J_active J(Y~P) + Molecules J_inactive J(Y~U) +end observables + +generate_network({overwrite=>1, max_stoich=>{J=>2}}); + +simulate_ode({t_end=>10000, n_steps=>10000,atol=>1e-08,rtol=>1e-08,sparse=>1}); diff --git a/tests/v1/bngl_corpus/Chattaraj_2021.bngl b/tests/v1/bngl_corpus/Chattaraj_2021.bngl new file mode 100644 index 00000000..5e53aefd --- /dev/null +++ b/tests/v1/bngl_corpus/Chattaraj_2021.bngl @@ -0,0 +1,74 @@ +# Clustering of three signaling molecules - Nephrin, Nck and NWASP +# References: +# 1. A. Chattaraj, M. L. Blinov and L. M. Loew. "The solubility product extends the buffering concept +# to heterotypic biomolecular condensates." eLife 2021, Vol. 10 Pages e67176, DOI: 10.7554/eLife.67176 (Figure 6) +# 2. A. Chattaraj and L. M. Loew. "The maximum solubility product marks the threshold for condensation +# of multivalent biomolecules". bioRxiv 2022, DOI: 10.1101/2022.10.04.510809 (Figure 6B) + +begin model + +begin parameters +# 1: Nephrin, 2: Nck, 3: NWASP +# Kd: binding affinity, kon: binding rate constant, koff: unbinding rate constant +kd_12 3500 +kd_23 3500 +koff_23 1000 +kon_23 koff_23/kd_23 +koff_12 1000 +kon_12 koff_12/kd_12 +end parameters + +begin molecule types +Nephrin(pY1,pY2,pY3) +Nck(S1,S2,S3,Sh2) +NWASP(p1,p2,p3,p4,p5,p6) +end molecule types + +begin seed species +1 Nephrin(pY1,pY2,pY3) 300 +2 Nck(S1,S2,S3,Sh2) 900 +3 NWASP(p1,p2,p3,p4,p5,p6) 450 +end seed species + +begin observables +Molecules tot_Nck Nck() +Molecules free_Nck Nck(S1,S2,S3,Sh2) +Molecules tot_NWASP NWASP() +Molecules free_NWASP NWASP(p1,p2,p3,p4,p5,p6) +Molecules tot_Nephrin Nephrin() +Molecules free_Nephrin Nephrin(pY1,pY2,pY3) +Molecules fully_bound_Nephrin Nephrin(pY1!+,pY2!+,pY3!+) +Molecules fully_bound_Nck Nck(S1!+,S2!+,S3!+,Sh2!+) +Molecules fully_bound_NWASP NWASP(p1!+,p2!+,p3!+,p4!+,p5!+,p6!+) +Molecules cluster_neph_nck_nw Nephrin().Nck().NWASP() +Molecules cluster_nck_nw Nck().NWASP() +end observables + + +begin reaction rules +_01_S1_P1: Nck(S1) + NWASP(p1) <-> Nck(S1!1).NWASP(p1!1) kon_23, koff_23 +_02_S2_P1: Nck(S2) + NWASP(p1) <-> Nck(S2!1).NWASP(p1!1) kon_23, koff_23 +_03_S3_P1: Nck(S3) + NWASP(p1) <-> Nck(S3!1).NWASP(p1!1) kon_23, koff_23 +_06_S3_P2: Nck(S3) + NWASP(p2) <-> Nck(S3!1).NWASP(p2!1) kon_23, koff_23 +_05_S2_P2: Nck(S2) + NWASP(p2) <-> Nck(S2!1).NWASP(p2!1) kon_23, koff_23 +_04_S1_P2: Nck(S1) + NWASP(p2) <-> Nck(S1!1).NWASP(p2!1) kon_23, koff_23 +_08_S2_P3: Nck(S2) + NWASP(p3) <-> Nck(S2!1).NWASP(p3!1) kon_23, koff_23 +_07_S1_P3: Nck(S1) + NWASP(p3) <-> Nck(S1!1).NWASP(p3!1) kon_23, koff_23 +_09_S3_P3: Nck(S3) + NWASP(p3) <-> Nck(S3!1).NWASP(p3!1) kon_23, koff_23 +_10_S1_P4: Nck(S1) + NWASP(p4) <-> Nck(S1!1).NWASP(p4!1) kon_23, koff_23 +_11_S2_P4: Nck(S2) + NWASP(p4) <-> Nck(S2!1).NWASP(p4!1) kon_23, koff_23 +_12_S3_P4: Nck(S3) + NWASP(p4) <-> Nck(S3!1).NWASP(p4!1) kon_23, koff_23 +_13_S1_P5: Nck(S1) + NWASP(p5) <-> Nck(S1!1).NWASP(p5!1) kon_23, koff_23 +_14_S2_P5: Nck(S2) + NWASP(p5) <-> Nck(S2!1).NWASP(p5!1) kon_23, koff_23 +_15_S3_P5: Nck(S3) + NWASP(p5) <-> Nck(S3!1).NWASP(p5!1) kon_23, koff_23 +_16_S1_P6: Nck(S1) + NWASP(p6) <-> Nck(S1!1).NWASP(p6!1) kon_23, koff_23 +_17_S3_P6: Nck(S3) + NWASP(p6) <-> Nck(S3!1).NWASP(p6!1) kon_23, koff_23 +_18_S2_P6: Nck(S2) + NWASP(p6) <-> Nck(S2!1).NWASP(p6!1) kon_23, koff_23 +_19_pY1_sh2: Nephrin(pY1) + Nck(Sh2) <-> Nephrin(pY1!1).Nck(Sh2!1) kon_12, koff_12 +_20_pY2_sh2: Nephrin(pY2) + Nck(Sh2) <-> Nephrin(pY2!1).Nck(Sh2!1) kon_12, koff_12 +_21_pY3_sh2: Nephrin(pY3) + Nck(Sh2) <-> Nephrin(pY3!1).Nck(Sh2!1) kon_12, koff_12 +end reaction rules + +end model + +writeXML() diff --git a/tests/v1/bngl_corpus/LR.bngl b/tests/v1/bngl_corpus/LR.bngl new file mode 100644 index 00000000..2ca8a6ee --- /dev/null +++ b/tests/v1/bngl_corpus/LR.bngl @@ -0,0 +1,43 @@ +## title: LR.bngl +## Simple ligand-receptor binding model expressed in cell units: +# time - seconds +# concentration - molecule number +# length - um, area - um^2, volume - um^3 +# bimolecular association - um^3/s (3D), um^2/s (2D) +## author: Jim Faeder +## date: 05Mar2018 + +begin model +begin parameters + NaV 6.02e8 # Conversion constant: M -> #/um^3 + Vcell 1000 # Typical eukaryotic cell volume ~ 1000 um^3 + Vec 1000*Vcell # Volume of extracellular space around each cell (1/cell density) + lig_conc 1e-8 # Ligand concentration - molar + L0 lig_conc*NaV*Vec # number of ligand molecules + R0 10000 # number of receptor molecules per cell + + kp1 1e6/(NaV*Vec) # Forward binding rate constant for L-R + km1 0.01 # Reverse binding rate constant for L-R +end parameters + +begin molecule types + L(r) # L molecule has one binding site for R + R(l) # R molecule has one binding site for L +end molecule types + +begin species + L(r) L0 + R(l) R0 +end species + +begin observables + Molecules FreeR R(l) + Molecules Bound L(r!1).R(l!1) +end observables + +begin reaction rules + L(r) + R(l) <-> L(r!1).R(l!1) kp1, km1 +end reaction rules +end model + +simulate({method=>"ode",t_end=>300,n_steps=>500}) diff --git a/tests/v1/bngl_corpus/LRR.bngl b/tests/v1/bngl_corpus/LRR.bngl new file mode 100644 index 00000000..1dc9aca9 --- /dev/null +++ b/tests/v1/bngl_corpus/LRR.bngl @@ -0,0 +1,53 @@ +## title: LRR.bngl +## Ligand-receptor binding and dimerization model expressed in cell units: +# time - seconds +# concentration - molecule number +# length - um, area - um^2, volume - um^3 +# bimolecular association - um^3/s (3D), um^2/s (2D) +## author: Jim Faeder +## date: 05Mar2018 + +begin model +begin parameters + NaV 6.02e8 # Conversion constant: M -> #/um^3 + Vcell 1000 # Typical eukaryotic cell volume ~ 1000 um^3 + Vec 1000*Vcell # Volume of extracellular space around each cell (1/cell density) + d_pm 0.01 # Effective thickness of the plasma membrane (10 nm) + Acell 1000 # Approximate area of PM + Vpm Acell*d_pm # Effective volume of PM + lig_conc 1e-8 # Ligand concentration - molar + L0 lig_conc*NaV*Vec # number of ligand molecules + R0 10000 # number of receptor molecules per cell + + kp1 1e6/NaV # #/um^3 1/s: Forward binding rate constant for L-R + km1 0.01 # 1/s Reverse binding rate constant for L-R +end parameters + +begin molecule types + L(r,r) # L molecule has two binding sites for R + R(l) # R molecule has one binding site for L +end molecule types + +begin species + L(r,r) L0/2 # L0 represents # of sites as opposed to molecules + R(l) R0 +end species + +begin observables + Molecules FreeR R(l) + Molecules Bound L(r,r!1).R(l!1) + Species Dimers R(l!1).L(r!1,r!2).R(l!2) +end observables + +begin reaction rules +LBind: L(r,r) + R(l) <-> L(r,r!1).R(l!1) kp1/Vec, km1 +Dimer: R(l) + L(r,r!1).R(l!1) <-> R(l!2).L(r!2,r!1).R(l!1) kp1/(Acell*d_pm), km1 +end reaction rules +end model + +simulate({method=>"ode",t_end=>300,n_steps=>500}) + +#parameter_scan({method=>"ode",parameter=>"lig_conc",par_min=>1e-12,par_max=>1e-6,\ +# n_scan_pts=>50,log_scale=>1,t_end=>1000,n_steps=>2}) + +#visualize({type=>"regulatory",background=>0,ruleNames=>1}) diff --git a/tests/v1/bngl_corpus/Motivating_example_cBNGL.bngl b/tests/v1/bngl_corpus/Motivating_example_cBNGL.bngl new file mode 100644 index 00000000..a5e2680e --- /dev/null +++ b/tests/v1/bngl_corpus/Motivating_example_cBNGL.bngl @@ -0,0 +1,225 @@ +# Signal transduction with receptor internalization and transcriptional reg. # +# cBNGL code: justin.s.hogg@gmail.com # +# conception: Leonard A. Harris, Justin S. Hogg, James R. Faeder # +# 6 June 2009 # +# # +# Demonstrates features of cBNGL in a biologically relevant scenario. # +# A motivating example for Winter Simulation Conference 2009 invited paper. # + +begin model +begin parameters + nEndo 5 # mean number of endosomes + + vol_EC 20.0 # volumes + vol_CP 4.0 + vol_NU 1.0 + vol_EN 0.1*nEndo + + sa_PM 0.4 # membrane surface areas + sa_NM 0.1 + sa_EM 0.01*nEndo + + eff_width 1.0 # effective surface width + + L0 1000 # initial species counts (extensive units: quantity, not concentration) + R0 200 + TF0 200 + DNA0 2 + Im0 40 + NP0 4 + + kp_LR 0.1 # kinetic parameters (2nd order reaction params in vol/time units) + km_LR 1.0 + kp_LL 0.1 + km_LL 1.0 + k_R_endo 1.0 + k_recycle 0.1 + k_R_transphos 1.0 + k_R_dephos 0.1 + kp_R_TF 0.1 + km_R_TF 0.1 + kp_R_TFp 0.1 + km_R_TFp 10.0 + k_TF_transphos 1.0 + k_TF_dephos 1.0 + kp_TF_TF 0.1 + km_TF_TF 1.0 + kp_TF_p1 0.1 + km_TF_p1 1.0 + k_transcribe 1.0 + k_translate 1.0 + k_mRNA_to_CP 1.0 + k_mRNA_deg 1.0 + k_P_deg 0.1 + k_Im_bind_CP 0.1 + k_Im_unbind_CP 0.1 + k_Im_bind_NU 0.1 + k_Im_unbind_NU 10.0 + k_Im_enters_NP 0.1 + k_Im_exits_NP 1.0 + k_Im_cross_NP 1.0 + kp_P1_p2 0.1 + km_P1_p2 1.0 +end parameters + +begin compartments + EC 3 vol_EC + PM 2 sa_PM * eff_width EC + CP 3 vol_CP PM + NM 2 sa_NM * eff_width CP + NU 3 vol_NU NM + EM 2 sa_EM * eff_width CP + EN 3 vol_EN EM +end compartments + +begin molecule types + L(r,d) # Ligand w/ receptor binding and dimerization sites. + R(l,tf~Y~pY) # Receptor with ligand and TF binding sites. + TF(r,d~Y~pY,dna,im) # Transcription factor (monomer) with receptor, DNA, and importin binding sites; and dimerization domain. + DNA(p1,p2) # DNA molecule with two promoter sites. + mRNA1() # mRNA transcript for Protein 1. + mRNA2() # mRNA transcript for Protein 2. + P1(im,dna) # Protein 1 with importin and DNA binding domains. + P2() # Protein 2. + Im(fg,cargo) # nuclear importin molecule with hydrophobic domain (fg) that interacts with nuclear pore. + NP(fg) # nuclear pore complex w/ hydrophobic FG repeat domain. + Sink() # a place for deleted molecules. +end molecule types + +begin species + L(r,d)@EC L0 + R(l,tf~Y)@PM R0 + TF(r,d~Y,dna,im)@CP TF0 + DNA(p1,p2)@NU DNA0 + Im(fg,cargo)@CP Im0 + NP(fg)@NM NP0 + + # Arbitrarily assign compartment CP to the abstract molecule "Sink". + $Sink()@CP 0 +end species + +begin reaction rules + # receptor-ligand binding. + Rule1: L(r) + R(l) <-> L(r!1).R(l!1) kp_LR, km_LR + + # ligand dimerization. + Rule2: L(d) + L(d) <-> L(d!1).L(d!1) kp_LL, km_LL + + # Rule3: receptor-dimer internalization. + Rule3: @PM:R.R -> @EM:R.R k_R_endo + + # receptor, ligand recycling. + Rule4: @EM:R -> @PM:R k_recycle + Rule5: @EN:L -> @EC:L k_recycle + + # receptor transphosphorylation. + Rule6: R.R(tf~Y) -> R.R(tf~pY) k_R_transphos + + # receptor dephosphorylation. + Rule7: R(tf~pY) -> R(tf~Y) k_R_dephos + + # receptor-TF binding. favor binding if TF(dim~Y), unbinding if TF(dim~pY). + Rule8: R(tf~pY) + TF(d~Y,r) <-> R(tf~pY!1).TF(d~Y,r!1) kp_R_TF, km_R_TF + Rule9: R(tf~pY) + TF(d~pY,r) <-> R(tf~pY!1).TF(d~pY,r!1) kp_R_TFp, km_R_TFp + + # transcription factor trans-phosphorylation. + Rule10: TF.R.R.TF(d~Y) -> TF.R.R.TF(d~pY) k_TF_transphos + + # transcription factor dephosphorylation (CP only). + Rule11: TF(d~pY)@CP -> TF(d~Y)@CP k_TF_dephos + + # transcription factor dimerization. + Rule12: TF(r,d~pY,dna) + TF(r,d~pY,dna) <-> TF(r,d~pY!1,dna).TF(r,d~pY!1,dna) kp_TF_TF, km_TF_TF + + # TF dimer binds promoter 1. + Rule13: TF(dna,im).TF(dna,im) + DNA(p1) <-> TF(dna!1,im).TF(dna!2,im).DNA(p1!1!2) kp_TF_p1, km_TF_p1 + + # transcription. + Rule14: DNA(p1!+) -> DNA(p1!+) + mRNA1()@NU k_transcribe + Rule15: DNA(p2!+) -> DNA(p2!+) + mRNA2()@NU k_transcribe + + # mRNA transport to cytoplams. + Rule16: mRNA1@NU -> mRNA1@CP k_mRNA_to_CP + Rule17: mRNA2@NU -> mRNA2@CP k_mRNA_to_CP + + # mRNA translation to protein. + Rule18: mRNA1@CP -> mRNA1@CP + P1(im,dna)@CP k_translate + Rule19: mRNA2@CP -> mRNA2@CP + P2()@CP k_translate + + # mRNA degradation. + Rule20: mRNA1 -> Sink()@CP k_mRNA_deg DeleteMolecules + Rule21: mRNA2 -> Sink()@CP k_mRNA_deg DeleteMolecules + + # protein degradation. + Rule22: P1 -> Sink()@CP k_P_deg DeleteMolecules + Rule23: P2 -> Sink()@CP k_P_deg DeleteMolecules + + # importin binds TF dimer (tends to pick up in CP, drop off in NU). + Rule24: TF(im,dna,r).TF(im,dna,r) + Im(cargo)@CP <-> TF(im!1,dna,r).TF(im!2,dna,r).Im(cargo!1!2)@CP k_Im_bind_CP, k_Im_unbind_CP + Rule25: TF(im,dna,r).TF(im,dna,r) + Im(cargo)@NU <-> TF(im!1,dna,r).TF(im!2,dna,r).Im(cargo!1!2)@NU k_Im_bind_NU, k_Im_unbind_NU + + # importin binds P1 (tends to pick up in CP, drop off in NU). + Rule26: P1(im,dna) + Im(cargo)@CP <-> P1(im!1,dna).Im(cargo!1)@CP k_Im_bind_CP, k_Im_unbind_CP + Rule27: P1(im,dna) + Im(cargo)@NU <-> P1(im!1,dna).Im(cargo!1)@NU k_Im_bind_NU, k_Im_unbind_NU + + # importin enters nuclear pore. + Rule28: Im(fg) + NP(fg) <-> Im(fg!1).NP(fg!1) k_Im_enters_NP, k_Im_exits_NP + + # importin traverses nuclear pore (with any cargo). + Rule29: Im(fg!1)@CP.NP(fg!1) <-> Im(fg!1)@NU.NP(fg!1) k_Im_cross_NP, k_Im_cross_NP MoveConnected + + # P1 binds promoter 2. + Rule30: P1(im,dna) + DNA(p2) <-> P1(im,dna!1).DNA(p2!1) kp_P1_p2, km_P1_p2 +end reaction rules + +begin observables + Molecules Tot_L L + Molecules Tot_R R + Molecules Tot_TF TF + Molecules Tot_DNA DNA + Molecules Tot_mRNA1 mRNA1 + Molecules Tot_mRNA2 mRNA2 + Molecules Tot_P1 P1 + Molecules Tot_P2 P2 + Molecules Tot_NP NP + Molecules Tot_Im Im + + Species L_Dimers_EC @EC:L.L + Species L_Dimers_PM @PM:L.L + Species L_Dimers_EN @EN:L.L + Species L_Dimers_EM @EM:L.L + + Molecules L_Bound_PM @PM:L + Molecules L_Bound_EM @EM:L + + Species R_Dimers_PM @PM:R.R + Species R_Dimers_EM @EM:R.R + + Molecules Catalytic_R R(tf~pY!?) + Molecules Catalytic_TF R(tf~pY!1).TF(r!1) + Molecules Phos_TF TF(d~pY!?) + + Species TF_Dimer_CP TF(d~pY!1)@CP.TF(d~pY!1)@CP + Species TF_Dimer_NU TF(d~pY!1)@NU.TF(d~pY!1)@NU + + Species Bound_prom1 DNA(p1!+) + Species Bound_prom2 DNA(p2!+) + + Species P1_NU P1@NU + Species P1_CP P1@CP + + Species Im_NU Im@NU + Species Im_CP Im@CP + + Species Im_Cargo_NP Im(fg!+,cargo!+) + + Species P1_NU_free P1(im,dna)@NU + Species P1_NU_dna P1(im,dna!+)@NU + + Species CountSink Sink()@CP +end observables +end model + +# actions # +generate_network({overwrite=>1}) +simulate({method=>"ode",t_start=>0,t_end=>40,n_steps=>100}) diff --git a/tests/v1/bngl_corpus/README.md b/tests/v1/bngl_corpus/README.md new file mode 100644 index 00000000..ae591d4f --- /dev/null +++ b/tests/v1/bngl_corpus/README.md @@ -0,0 +1,36 @@ +# Curated BNGL corpus fixtures + +Public community BNGL models, chosen for feature coverage, backing +`tests/v1/test_bngl_corpus.py`: it asserts `parse_bngl` enumerates the same +model entities BNG2.pl does. The BNG2.pl answers are cached in `golden.json` +(the entity name sets BNG2.pl emits from `writeModel`, its canonical parse), +so the test needs **no BNG2.pl** -- it compares the reader against the frozen +oracle. Regenerate the golden (needs BNG2.pl) with: + + python tests/v1/test_bngl_corpus.py + +Vendored from public repos (RuleWorld/RuleHub, wshlavacek/BNGL-Models). + +| model | source | upstream path | exercises | +|---|---|---|---| +| `An_2009.bngl` | rulehub | `Published/An2009/An_2009.bngl` | labeled seed species; large rule-based (TLR signaling) | +| `Barua_2009__PATCHED.bngl` | rulehub | `Published/Barua2009/Barua_2009.bngl` | uses the `begin species` block alias | +| `Chattaraj_2021.bngl` | rulehub | `Published/Chattaraj2021/Chattaraj_2021.bngl` | indexed seed species | +| `LR.bngl` | rulehub | `Tutorials/NativeTutorials/LR/LR.bngl` | uses the `begin species` block alias | +| `LRR.bngl` | rulehub | `Tutorials/NativeTutorials/LRR/LRR.bngl` | uses the `begin species` block alias | +| `Motivating_example_cBNGL.bngl` | rulehub | `Tutorials/MotivatingexamplecBNGL/Motivating_example_cBNGL.bngl` | compartmental BNGL | +| `Ras_bistability_v2.bngl` | bngl_models | `my_models/ode/Ras_bistability_v2.bngl` | component reordering vs BNG canonical order | +| `Rule_based_egfr_tutorial.bngl` | rulehub | `Published/Rulebasedegfrtutorial/Rule_based_egfr_tutorial.bngl` | indexed seed species | +| `akt-signaling.bngl` | rulehub | `Examples/biology/aktsignaling/akt-signaling.bngl` | states + bonds; multisite phosphorylation | +| `apoptosis-cascade.bngl` | rulehub | `Examples/biology/apoptosiscascade/apoptosis-cascade.bngl` | states + bonds; compartment-free signaling | +| `bcr-signaling.bngl` | rulehub | `Examples/biology/bcrsignaling/bcr-signaling.bngl` | $ clamp; states + bonds (B-cell receptor) | +| `blood-coagulation-thrombin.bngl` | rulehub | `Examples/biology/bloodcoagulationthrombin/blood-coagulation-thrombin.bngl` | large rule-based; local functions | +| `bmp-signaling.bngl` | rulehub | `Examples/biology/bmpsignaling/bmp-signaling.bngl` | complex formation; $ clamp | +| `brusselator-oscillator.bngl` | rulehub | `Examples/biology/brusselatoroscillator/brusselator-oscillator.bngl` | small sanity model | +| `catalysis.bngl` | bngl_models | `my_models/ode/catalysis.bngl` | compartmental BNGL; $ clamp under @compartment prefix | +| `egg.bngl` | rulehub | `Published/Hlavacek2018Egg/egg.bngl` | line continuations; bare-molecule seed (t -> t()) | +| `elephant_EFA.bngl` | rulehub | `Published/Hlavacek2018Elephant/elephant_EFA.bngl` | bare-molecule seed; continuations | +| `energy_transport_pump.bngl` | rulehub | `Examples/energy/energytransportpump/energy_transport_pump.bngl` | energy patterns | +| `example1.bngl` | rulehub | `Tutorials/example1/example1.bngl` | indexed parameters + indexed seed species | +| `genetic_bistability_energy.bngl` | rulehub | `Examples/genetics/geneticbistabilityenergy/genetic_bistability_energy.bngl` | energy patterns | +| `immob_equiv_lig_sites.bngl` | bngl_models | `my_models/nf/immob_equiv_lig_sites.bngl` | line continuations in functions; local functions | diff --git a/tests/v1/bngl_corpus/Ras_bistability_v2.bngl b/tests/v1/bngl_corpus/Ras_bistability_v2.bngl new file mode 100644 index 00000000..84f47c6c --- /dev/null +++ b/tests/v1/bngl_corpus/Ras_bistability_v2.bngl @@ -0,0 +1,421 @@ +# A model for Ras signaling +# This model accounts for two compartments: +# 1) plasma membrane and 2) cytosol +# +# Author information: +# William S. Hlavacek +# Theoretical Division +# [PII REMOVED] +# +# Date: +# 26-December-2013 +# +# Platform information: +# Mac OS X, version 10.5.8 (desktop), version 10.6.8 (laptop) +# +# Software information: +# BioNetGen, version 2.2.5 (http://bionetgen.org) +# NFsim, version 1.11 + + +begin model + + +begin parameters + +# simulation parameters + +# fraction of a single cell to be considered in a stochastic simulation +f 1.0 # [=] dimensionless (0\ +RTK(pTyrGEF!1).SOS1(GEF,REM,GRB2_SH2!1) kaGEF_cyto +# RTK binds plasma membrane-tethered SOS1 +# SOS1 is tethered at GEF only +RTK(pTyrGEF)+SOS1(GEF!+,REM,GRB2_SH2)->\ +RTK(pTyrGEF!1).SOS1(GEF!+,REM,GRB2_SH2!1) kaGEF_pm +# SOS1 is tethered at REM only +RTK(pTyrGEF)+SOS1(GEF,REM!+,GRB2_SH2)->\ +RTK(pTyrGEF!1).SOS1(GEF,REM!+,GRB2_SH2!1) kaGEF_pm +# SOS1 is tethered at both GEF and REM +RTK(pTyrGEF)+SOS1(GEF!+,REM!+,GRB2_SH2)->\ +RTK(pTyrGEF!1).SOS1(GEF!+,REM!+,GRB2_SH2!1) kaGEF_pm +# reverse reaction +RTK(pTyrGEF!1).SOS1(GRB2_SH2!1)->RTK(pTyrGEF)+SOS1(GRB2_SH2) kdGEF + +# reversible binding of GAP to RTK +# forward reaction +# RTK binds cytosolic RASA1 +RTK(pTyrGAP)+RASA1(tSH2,GAP)->RTK(pTyrGAP!1).RASA1(tSH2!1,GAP) kaGAP_cyto +# RTK binds plasma membrane-tethered RASA1 +RTK(pTyrGAP)+RASA1(tSH2,GAP!+)->RTK(pTyrGAP!1).RASA1(tSH2!1,GAP!+) kaGAP_pm +# reverse reaction +RTK(pTyrGAP!1).RASA1(tSH2!1)->RTK(pTyrGAP)+RASA1(tSH2) kdGAP + +# To do: add rules for REM effect +# To do: add rules for KRASG12D + +# GDP binds the apo form of Ras +KRAS(G~apo)+GDP->KRAS(G~GDP) kp1 +# GDP spontaneously dissociates from Ras +KRAS(G~GDP)->KRAS(G~apo)+GDP km1 + +# GEF reversibly binds apo Ras +# cytosolic GEF binds apo Ras +KRAS(G~apo)+SOS1(GEF,REM,GRB2_SH2)->\ +KRAS(G~apo!1).SOS1(GEF!1,REM,GRB2_SH2) kp2_cyto +# plasma membrane-recruited GEF binds apo Ras +# SOS1 is tethered at REM only +KRAS(G~apo)+SOS1(GEF,REM!+,GRB2_SH2)->\ +KRAS(G~apo!1).SOS1(GEF!1,REM!+,GRB2_SH2) kp2_pm +# SOS1 is tethered at GRB2_SH2 only +KRAS(G~apo)+SOS1(GEF,REM,GRB2_SH2!+)->\ +KRAS(G~apo!1).SOS1(GEF!1,REM,GRB2_SH2!+) kp2_pm +# SOS1 is tethered at both REM and GRB2_SH2 +KRAS(G~apo)+SOS1(GEF,REM!+,GRB2_SH2!+)->\ +KRAS(G~apo!1).SOS1(GEF!1,REM!+,GRB2_SH2!+) kp2_pm +# reverse reaction +KRAS(G~apo!1).SOS1(GEF!1)->KRAS(G~apo)+SOS1(GEF) km2 + +# GEF reversibly binds RasGDP +# cytosolic GEF binds RasGDP +KRAS(G~GDP)+SOS1(GEF,REM,GRB2_SH2)->\ +KRAS(G~GDP!1).SOS1(GEF!1,REM,GRB2_SH2) kp3_cyto +# plasma membrane-recruited GEF binds RasGDP +# SOS1 is tethered at REM only +KRAS(G~GDP)+SOS1(GEF,REM!+,GRB2_SH2)->\ +KRAS(G~GDP!1).SOS1(GEF!1,REM!+,GRB2_SH2) kp3_pm +# SOS1 is tethered at GRB2_SH2 only +KRAS(G~GDP)+SOS1(GEF,REM,GRB2_SH2!+)->\ +KRAS(G~GDP!1).SOS1(GEF!1,REM,GRB2_SH2!+) kp3_pm +# SOS1 is tethered at both REM and GRB2_SH2 +KRAS(G~GDP)+SOS1(GEF,REM!+,GRB2_SH2!+)->\ +KRAS(G~GDP!1).SOS1(GEF!1,REM!+,GRB2_SH2!+) kp3_pm +# reverse reaction +KRAS(G~GDP!1).SOS1(GEF!1)->KRAS(G~GDP)+SOS1(GEF) km3 + +# GDP binds apo Ras-GEF complex +KRAS(G~apo!1).SOS1(GEF!1)+GDP->KRAS(G~GDP!1).SOS1(GEF!1) kp4 +# GDP dissociates from RasGDP-GEF complex +KRAS(G~GDP!1).SOS1(GEF!1)->KRAS(G~apo!1).SOS1(GEF!1)+GDP km4 + +# GTP binds the apo form of Ras +KRAS(G~apo)+GTP->KRAS(G~GTP) kp5 +# GTP spontaneoulsy dissociates from Ras +KRAS(G~GTP)->KRAS(G~apo)+GTP km5 + +# GTP binds apo Ras-GEF complex +KRAS(G~apo!1).SOS1(GEF!1)+GTP->KRAS(G~GTP!1).SOS1(GEF!1) kp6 +# GTP dissociates from RasGTP-GEF complex +KRAS(G~GTP!1).SOS1(GEF!1)->KRAS(G~apo!1).SOS1(GEF!1)+GTP km6 + +# GEF reversibly binds RasGTP +# cytosolic GEF binds RasGTP +KRAS(G~GTP)+SOS1(GEF,REM,GRB2_SH2)->\ +KRAS(G~GTP!1).SOS1(GEF!1,REM,GRB2_SH2) kp7_cyto +# plasma membrane-recruited GEF binds RasGTP +# SOS1 is tethered at REM only +KRAS(G~GTP)+SOS1(GEF,REM!+,GRB2_SH2)->\ +KRAS(G~GTP!1).SOS1(GEF!1,REM!+,GRB2_SH2) kp7_pm +# SOS1 is tethered at GRB2_SH2 only +KRAS(G~GTP)+SOS1(GEF,REM,GRB2_SH2!+)->\ +KRAS(G~GTP!1).SOS1(GEF!1,REM,GRB2_SH2!+) kp7_pm +# SOS1 is tethered at both REM and GRB2_SH2 +KRAS(G~GTP)+SOS1(GEF,REM!+,GRB2_SH2!+)->\ +KRAS(G~GTP!1).SOS1(GEF!1,REM!+,GRB2_SH2!+) kp7_pm +# reverse reaction +KRAS(G~GTP!1).SOS1(GEF!1)->KRAS(G~GTP)+SOS1(GEF) km7 + +# Effector reversibly binds RasGTP +KRAS(G~GTP)+CRAF(RBD)<->KRAS(G~GTP!1).CRAF(RBD!1) kaEff,kdEff + +# GTP hydrolysis via the intrinsic GTPase activity of Ras +# RasGTP is free +KRAS(G~GTP)->KRAS(G~GDP) khyd +# RasGTP is bound to effector +# fast release of effector after hydrolysis +KRAS(G~GTP!1).CRAF(RBD!1)->KRAS(G~GDP)+CRAF(RBD) khyd + +# GAP-assisted hydrolysis of GTP via a Michaelis-Menten mechanism +# Enzyme-substrate association +# GAP is cytosolic +KRAS(G~GTP)+RASA1(tSH2,GAP)->KRAS(G~GTP!1).RASA1(tSH2,GAP!1) kf_cyto +# GAP is plasma membrane recruited +KRAS(G~GTP)+RASA1(tSH2!+,GAP)->KRAS(G~GTP!1).RASA1(tSH2!+,GAP!1) kf_pm +# Enzyme-substrate dissociation (no reaction) +KRAS(G~GTP!1).RASA1(GAP!1)->KRAS(G~GTP)+RASA1(GAP) kr +# Hydrolysis (reaction) +# fast dissoc of enzyme and product after hydrolysis +KRAS(G~GTP!1).RASA1(GAP!1)->KRAS(G~GDP)+RASA1(GAP) kcat + +end reaction rules + + +end model + + +# actions + +generate_network({overwrite=>1}); + +# equilibration +# simulate 6 days, report progress every 4 hours +simulate_ode({suffix=>"equil",t_start=>-518400,t_end=>0,n_steps=>36}) +saveConcentrations() + +# simulate response to GEF & GAP recruitment to plasma membrane + +resetConcentrations() +setConcentration("RTK(pTyrGEF,pTyrGAP)",1.0e3) +simulate_ode({suffix=>"1",t_start=>0,t_end=>120,n_steps=>120}) + +resetConcentrations() +setConcentration("RTK(pTyrGEF,pTyrGAP)",5.0e3) +simulate_ode({suffix=>"2",t_start=>0,t_end=>120,n_steps=>120}) + +resetConcentrations() +setConcentration("RTK(pTyrGEF,pTyrGAP)",1.0e4) +simulate_ode({suffix=>"3",t_start=>0,t_end=>120,n_steps=>120}) + +resetConcentrations() +setConcentration("RTK(pTyrGEF,pTyrGAP)",5.0e4) +simulate_ode({suffix=>"4",t_start=>0,t_end=>120,n_steps=>120}) + +resetConcentrations() +setConcentration("RTK(pTyrGEF,pTyrGAP)",1.0e5) +simulate_ode({suffix=>"5",t_start=>0,t_end=>120,n_steps=>120}) diff --git a/tests/v1/bngl_corpus/Rule_based_egfr_tutorial.bngl b/tests/v1/bngl_corpus/Rule_based_egfr_tutorial.bngl new file mode 100644 index 00000000..cce94a8a --- /dev/null +++ b/tests/v1/bngl_corpus/Rule_based_egfr_tutorial.bngl @@ -0,0 +1,55 @@ +begin model + +begin compartments +c0 3 1 +end compartments + +begin parameters +end parameters + +begin molecule types +EGF(Site) +EGFR(ecd,tmd,Y1~u~p,Y2~u~p) +Grb2(sh2) +Shc(sh3,Y~u~p) +end molecule types + +begin seed species +1 @c0:EGFR(ecd,tmd,Y1~u,Y2~u) 100.0 +2 @c0:EGF(Site) 680.0 +3 @c0:Grb2(sh2) 58.0 +4 @c0:Shc(sh3,Y~p) 0.0 +5 @c0:Shc(sh3,Y~u) 150.0 +end seed species + +begin observables +Molecules O0_EGF_tot @c0:EGF() +Molecules O0_EGFR_tot @c0:EGFR() +Molecules O0_Grb2_tot @c0:Grb2() +Molecules O0_Shc_tot @c0:Shc() +Molecules Dimers @c0:EGFR(tmd!+) +Species Dimers_s @c0:EGFR(tmd!+) +Molecules Y1 @c0:EGFR(Y1~p!?) +Molecules Y2 @c0:EGFR(Y2~p!?) +Molecules Y_total @c0:EGFR(Y1~p!?) @c0:EGFR(Y2~p!?) +end observables + +begin functions +end functions + +begin reaction rules +ligand_bind: @c0:EGFR(ecd,tmd) + @c0:EGF(Site) <-> @c0:EGFR(ecd!1,tmd).EGF(Site!1) 0.003, 0.06 +dimeriz: @c0:EGFR(ecd!+,tmd)%1 + @c0:EGFR(ecd!+,tmd)%2 <-> @c0:EGFR(ecd!+,tmd!1)%1.EGFR(ecd!+,tmd!1)%2 0.001, 0.1 +Y2_phosph: @c0:EGFR(tmd!+,Y2~u) -> @c0:EGFR(tmd!+,Y2~p) 0.5 +Y1_phosph: @c0:EGFR(tmd!+,Y1~u) -> @c0:EGFR(tmd!+,Y1~p) 0.5 +Y2_dephosph: @c0:EGFR(Y2~p) -> @c0:EGFR(Y2~u) 4.5 +Y1_dephosph: @c0:EGFR(Y1~p) -> @c0:EGFR(Y1~u) 4.5 +R_Grb2_interaction: @c0:EGFR(Y1~p) + @c0:Grb2(sh2) <-> @c0:EGFR(Y1~p!1).Grb2(sh2!1) 0.001, 0.05 +R_ShcU_interaction: @c0:EGFR(Y2~p) + @c0:Shc(sh3,Y~u) <-> @c0:EGFR(Y2~p!1).Shc(sh3!1,Y~u) 0.045, 0.6 +Shc_phosph: @c0:EGFR(Y2~p!1).Shc(sh3!1,Y~u) -> @c0:EGFR(Y2~p!1).Shc(sh3!1,Y~p) 3.0 +R_ShcP_interaction: @c0:EGFR(Y2~p) + @c0:Shc(sh3,Y~p) <-> @c0:EGFR(Y2~p!1).Shc(sh3!1,Y~p) 4.5E-4, 0.3 +end reaction rules + +end model + +generate_network({max_iter=>12,max_agg=>12,max_stoich=>{EGF=>100,EGFR=>100,Grb2=>100,Shc=>100},overwrite=>1}) diff --git a/tests/v1/bngl_corpus/akt-signaling.bngl b/tests/v1/bngl_corpus/akt-signaling.bngl new file mode 100644 index 00000000..3945de17 --- /dev/null +++ b/tests/v1/bngl_corpus/akt-signaling.bngl @@ -0,0 +1,100 @@ +begin model +begin parameters + # Signaling rates + k_bind 1e-3 # RTK-GF binding + k_unbind 1e-4 # GF dissociation + k_rtk_act 0.5 # Autophosphorylation + k_pi3k 0.2 # PI3K recruitment/activation + k_akt_t308 0.5 # PDK1-mediated T308 phos + k_akt_s473 0.3 # mTORC2-mediated S473 phos + k_mtorc1 0.4 # Akt-mediated mTORC1 activation + k_s6k_act 0.2 # mTORC1-mediated S6K phos + + # Negative feedback / resetting + k_reset 0.1 # General phosphatase/deactivation + k_fb 0.05 # S6K-mediated inhibition of RTK signaling + + # Total concentrations + GF_tot 100 + RTK_tot 50 + PI3K_tot 40 + AKT_tot 100 + mTORC2_tot 20 + mTORC1_tot 50 + S6K_tot 30 +end parameters + +begin molecule types + GrowthFactor(r) + RTK(l,state~U~P,fb_site~open~closed) + PI3K(state~off~on) + # Akt has two critical phosphorylation sites + AKT(t308~U~P,s473~U~P) + mTORC2(state~on) + mTORC1(state~off~on) + S6K(state~U~P) +end molecule types + +begin seed species + GrowthFactor(r) GF_tot + RTK(l,state~U,fb_site~open) RTK_tot + PI3K(state~off) PI3K_tot + AKT(t308~U,s473~U) AKT_tot + mTORC2(state~on) mTORC2_tot + mTORC1(state~off) mTORC1_tot + S6K(state~U) S6K_tot +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules Active_RTK RTK(state~P) # Receptor tyrosine kinase activation + Molecules Active_PI3K PI3K(state~on) # PI3K recruitment to membrane + Molecules pAkt_T308 AKT(t308~P) # PDK1-mediated phosphorylation + Molecules pAkt_S473 AKT(s473~P) # mTORC2-mediated phosphorylation + Molecules Double_pAkt AKT(t308~P,s473~P) # Fully active Akt + Molecules Active_mTORC1 mTORC1(state~on) # Downstream growth effector +end observables + +begin reaction rules + ## RECEPTOR ACTIVATION & FEEDBACK + # Growth factor binds to open RTK + GrowthFactor(r) + RTK(l,fb_site~open) <-> GrowthFactor(r!1).RTK(l!1,fb_site~open) k_bind,k_unbind + + # Ligand-bound RTK becomes active + GrowthFactor(r!1).RTK(l!1,state~U) -> GrowthFactor(r!1).RTK(l!1,state~P) k_rtk_act + + ## PI3K / AKT PATHWAY + # Active RTK recruits and activates PI3K + RTK(state~P) + PI3K(state~off) -> RTK(state~P) + PI3K(state~on) k_pi3k + + # Active PI3K facilitates Akt phosphorylation at T308 (PDK1 proxy) + PI3K(state~on) + AKT(t308~U) -> PI3K(state~on) + AKT(t308~P) k_akt_t308 + + # Constitutive mTORC2 phosphorylates Akt at S473 + mTORC2() + AKT(s473~U) -> mTORC2() + AKT(s473~P) k_akt_s473 + + ## MTOR SIGNALING & FEEDBACK + # Double-phosphorylated Akt activates mTORC1 + AKT(t308~P,s473~P) + mTORC1(state~off) -> AKT(t308~P,s473~P) + mTORC1(state~on) k_mtorc1 + + # mTORC1 activates S6K,which mediates feedback + mTORC1(state~on) + S6K(state~U) -> mTORC1(state~on) + S6K(state~P) k_s6k_act + + # Negative feedback: pS6K inhibits RTK ligand binding by "closing" the feedback site + S6K(state~P) + RTK(fb_site~open) -> S6K(state~P) + RTK(fb_site~closed) k_fb + + ## RESETTING MECHANISMS + RTK(state~P) -> RTK(state~U) k_reset + PI3K(state~on) -> PI3K(state~off) k_reset + AKT(t308~P) -> AKT(t308~U) k_reset + AKT(s473~P) -> AKT(s473~U) k_reset + mTORC1(state~on) -> mTORC1(state~off) k_reset + S6K(state~P) -> S6K(state~U) k_reset + RTK(fb_site~closed) -> RTK(fb_site~open) k_reset +end reaction rules + +begin actions + generate_network({overwrite=>1}) + simulate({method=>"ode",t_start=>0,t_end=>500,n_steps=>250}) +end actions +end model diff --git a/tests/v1/bngl_corpus/apoptosis-cascade.bngl b/tests/v1/bngl_corpus/apoptosis-cascade.bngl new file mode 100644 index 00000000..6c3a00a1 --- /dev/null +++ b/tests/v1/bngl_corpus/apoptosis-cascade.bngl @@ -0,0 +1,110 @@ +begin model +begin parameters + # Apoptosis cascade: Integrated extrinsic and intrinsic death signaling. + # Advanced features: Hill kinetics for MOMP and MoveConnected for CytoC. + + # Upstream Induction + k_ligand 0.8 # Death ligand drive + + # Mitochondrial Decision (MOMP) + k_bid_trunc 1.0 # BID -> tBID + k_momp_max 5.0 # Final mitochondrial commitment + Km_tbid 150 # Threshold for MOMP + n_hill 4.0 # Sharp switching behavior + + # Execution + k_apaf_act 1.5 # CytoC/Apaf-1 assembly + k_c3_act 2.5 # Caspase-3 execution + + # Regulation (Survival) + k_xiap_bind 5.0 # XIAP sequesters C3 + k_smac_rec 2.0 # SMAC inhibits XIAP + + # Initials + C8_tot 500 + Bid_tot 400 + MOMP_sites 100 # Mitochondrial surface proxy + Apaf1_tot 300 + C3_tot 1000 + XIAP_tot 200 + SMAC_tot 150 +end parameters + +begin molecule types + DeathLigand() + Caspase8(s~U~A) + Bid(s~U~T) + Mito(state~intact~leaky,loc~mit~cyt) + Apaf1(s~U~A) + Caspase3(b,s~U~A) + XIAP(b1,b2) + SMAC(b,loc~mit~cyt) +end molecule types + +begin seed species + DeathLigand() 10 + Caspase8(s~U) C8_tot + Bid(s~U) Bid_tot + Mito(state~intact,loc~mit) MOMP_sites + Apaf1(s~U) Apaf1_tot + Caspase3(b,s~U) C3_tot + XIAP(b1,b2) XIAP_tot + SMAC(b,loc~mit) SMAC_tot +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules Global_Death Caspase3(s~A) # Executioner activity + Molecules Mitochondrial_Fail Mito(state~leaky) # Irreversible commitment + Molecules CytoC_Released Mito(loc~cyt) # Spatial readout + Molecules Inhibitor_Sequest Caspase3(b!1).XIAP(b1!1) # Survival status + Molecules Mitochondria_Gate Mito(state~leaky) # Commitment marker + Molecules Active_Bid Bid(s~T) # Feedback driver status +end observables + +begin functions + # Sharp switch for mitochondrial collapse (MOMP) + v_momp() = k_momp_max * (Active_Bid^n_hill) / (Km_tbid^n_hill + Active_Bid^n_hill) +end functions + +begin reaction rules + ## RECEPTION & INITIATION (Extrinsic) + # Death ligand triggers Caspase-8 + DeathLigand() + Caspase8(s~U) -> DeathLigand() + Caspase8(s~A) k_ligand + + # Caspase-8 truncates Bid + Caspase8(s~A) + Bid(s~U) -> Caspase8(s~A) + Bid(s~T) k_bid_trunc + + ## MITOCHONDRIAL COLLAPSE (Intrinsic) + # tBid triggers MOMP (Functional Rate - Hill) + Mito(state~intact) -> Mito(state~leaky) v_momp() + + # MoveConnected: CytoC and SMAC move to cytosol upon MOMP + Mito(state~leaky,loc~mit) -> Mito(state~leaky,loc~cyt) 10.0 MoveConnected + SMAC(loc~mit) + Mito(state~leaky) -> SMAC(loc~cyt) + Mito(state~leaky) 10.0 + + ## EXECUTION + # CytoC (proxied by leaky mito) activates Apaf-1 + Mito(state~leaky,loc~cyt) + Apaf1(s~U) -> Mito(state~leaky,loc~cyt) + Apaf1(s~A) k_apaf_act + + # Initiators (C8,Apaf-C9) activate Caspase-3 + Caspase8(s~A) + Caspase3(s~U) -> Caspase8(s~A) + Caspase3(s~A) k_c3_act + Apaf1(s~A) + Caspase3(s~U) -> Apaf1(s~A) + Caspase3(s~A) k_c3_act + + ## REGULATION & SURVIVAL + # XIAP inhibits active C3 + Caspase3(s~A,b) + XIAP(b1) <-> Caspase3(s~A,b!1).XIAP(b1!1) k_xiap_bind,0.1 + + # SMAC inhibits XIAP,letting C3 free + SMAC(loc~cyt,b) + XIAP(b2) <-> SMAC(loc~cyt,b!1).XIAP(b2!1) k_smac_rec,0.1 + + ## RESET (Baseline recovery - failed apoptosis case) + Caspase3(s~A) -> Caspase3(s~U) 0.05 + Apaf1(s~A) -> Apaf1(s~U) 0.05 +end reaction rules + +begin actions + generate_network({overwrite=>1}) + simulate({method=>"ode",t_end=>10,n_steps=>200}) +end actions +end model diff --git a/tests/v1/bngl_corpus/bcr-signaling.bngl b/tests/v1/bngl_corpus/bcr-signaling.bngl new file mode 100644 index 00000000..1fb5bfba --- /dev/null +++ b/tests/v1/bngl_corpus/bcr-signaling.bngl @@ -0,0 +1,111 @@ +begin model +begin parameters + # BCR signaling: The B-cell antigen receptor cascade. + # Advanced features: Hill kinetics for calcium and DeleteMolecules for clearance. + + # Recognition & Initiation + k_bind_ag 1e-4 # BCR-Antigen binding + k_lyn_phos 2.0 # Lyn-mediated ITAM phos + + # Signal Relay (Syk/PLCg2) + k_syk_act 1.5 # Syk recruitment and activation + k_plcg_act 1.0 # PLCgamma2 activation + + # Downstream Messengers (Calcium Phase) + k_ca_max 10.0 # Peak calcium flux + Km_ca 300 # Threshold for calcium response + n_hill 3.0 # Sharp trigger mechanics + + # Negative Feedback (CD22/SHP-1) + k_shp_rec 1.2 # SHP-1 recruitment to CD22 + k_dephos_bcr 5.0 # SHP-1 mediated dephosphorylation (Fast) + + # Clearance + k_endo 0.2 # Receptor internalisation + + # Initials + BCR_tot 500 + Antigen_tot 50 + Syk_tot 120 + PLCg2_tot 100 + CD22_tot 150 + SHP1_tot 150 + Ca_cyt 10 # Basal calcium +end parameters + +begin molecule types + BCR(Y~U~P,b) + Antigen(b) + Syk(s~U~A) + PLCg2(s~U~A) + CD22(Y~U~P,b) + SHP1(b) + Calcium() +end molecule types + +begin seed species + BCR(Y~U,b) BCR_tot + Antigen(b) Antigen_tot + Syk(s~U) Syk_tot + PLCg2(s~U) PLCg2_tot + CD22(Y~U,b) CD22_tot + SHP1(b) SHP1_tot + Calcium() Ca_cyt +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules Calcium_Signal Calcium() # Secondary relay status + Molecules Active_PLCg2 PLCg2(s~A) # Flux driver status + Molecules Active_Relay Syk(s~A) # Initiator kinase status + Molecules Effector_Load PLCg2(s~A) # Metabolic trigger level + Molecules Negative_Brake SHP1(b!+) # Feedback intensity + Molecules Sync_Antigen Antigen(b!+) # Engagement metrics +end observables + +begin functions + # Sharp calcium flux (Hill kinetics) + v_ca() = k_ca_max * (Active_PLCg2^n_hill) / (Km_ca^n_hill + Active_PLCg2^n_hill) +end functions + +begin reaction rules + ## ACTIVATION PHASE + # Antigen binds BCR and triggers ITAM phos (Syk platform) + BCR(b) + Antigen(b) <-> BCR(b!1).Antigen(b!1) k_bind_ag,0.1 + BCR(Y~U,b!+) -> BCR(Y~P,b!+) k_lyn_phos + + # p-BCR recruits and activates Syk + BCR(Y~P) + Syk(s~U) -> BCR(Y~P) + Syk(s~A) k_syk_act + + ## EFFECTOR PHASE + # Syk activates PLCgamma2 + Syk(s~A) + PLCg2(s~U) -> Syk(s~A) + PLCg2(s~A) k_plcg_act + + # PLCg2 triggers Calcium flux (Functional Rate) + 0 -> Calcium() v_ca() + + ## NEGATIVE FEEDBACK + # CD22 gets phosphorylated (Slow) + CD22(Y~U) -> CD22(Y~P) 0.1 + + # p-CD22 recruits SHP-1 + CD22(Y~P) + SHP1(b) <-> CD22(Y~P!1).SHP1(b!1) k_shp_rec,0.2 + + # SHP-1 dephosphorylates BCR ITAM (Brake) + SHP1(b!+) + BCR(Y~P) -> SHP1(b!+) + BCR(Y~U) k_dephos_bcr + + ## CLEARANCE + # DeleteMolecules: Clear signaled receptors + BCR(b!1).Antigen(b!1) -> 0 k_endo DeleteMolecules + + ## RESET + Calcium() -> 0 0.2 # Calcium pumps + Syk(s~A) -> Syk(s~U) 0.1 + PLCg2(s~A) -> PLCg2(s~U) 0.1 +end reaction rules + +begin actions + generate_network({overwrite=>1}) + simulate({method=>"ode",t_end=>80,n_steps=>160}) +end actions +end model diff --git a/tests/v1/bngl_corpus/blood-coagulation-thrombin.bngl b/tests/v1/bngl_corpus/blood-coagulation-thrombin.bngl new file mode 100644 index 00000000..ee432c50 --- /dev/null +++ b/tests/v1/bngl_corpus/blood-coagulation-thrombin.bngl @@ -0,0 +1,101 @@ +begin model +begin parameters + # Blood coagulation: Thrombin burst and feedback propagation. + # Advanced features: TotalRate for clotting and DeleteMolecules for inhibition. + + # Initiation (Tissue Factor Pathway) + k_tf_vii 0.1 # Tissue Factor + FVIIa + k_xa_act 1.0 # Initial Xa generation + + # Amplification Loop (Thrombin Burst) + k_v_act 1.2 # Thrombin activates FV + k_viii_act 1.0 # Thrombin activates FVIII + k_complex_va_xa 5.0 # Prothrombinase assembly + k_pt_burst 50.0 # Efficient Thrombin generation by complex + + # Fibrin Formation + k_fibrin_synth 2.0 # Thrombin converts Fibrinogen + Km_clot 500 # Surface limit + + # Anticoagulation (Neutralization) + k_at_inh 1.5 # Antithrombin-III neutralizing factors + k_tpi_inh 1.0 # TFPI inhibition + k_deg_complex 10.0 # Rapid removal of inhibited complexes + + # Initials + Prothrom_tot 1200 + Fibrino_tot 2000 + FactorX_tot 100 + FactorV_tot 50 + ATIII_tot 300 + TF_sites 20 # Vascular injury signal +end parameters + +begin molecule types + TF(b,s~active) + FactorX(b,s~U~A) + FactorV(b,s~U~A) + Prothrombin(s~U~A) + Thrombin(s~U~A) + Fibrinogen(s~S~F) # Soluble vs Fibrin + AT(b) # Antithrombin +end molecule types + +begin seed species + TF(b,s~active) TF_sites + FactorX(b,s~U) FactorX_tot + FactorV(b,s~U) FactorV_tot + Prothrombin(s~U) Prothrom_tot + Fibrinogen(s~S) Fibrino_tot + AT(b) ATIII_tot +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules Clot_Mass Fibrinogen(s~F) # Final functional output + Molecules Active_Thrombin Thrombin(s~A) # Feedback driver + Molecules Proth_Complex FactorX(b!1).FactorV(b!1) # Assembly platform + Molecules Feedback_Factors FactorV(s~A) # Priming status + Molecules Inhibitor_Work AT(b!+) # Clearance efficiency +end observables + +begin reaction rules + ## INITIATION Phase + # Tissue Factor activates Factor X + TF(s~active) + FactorX(s~U) -> TF(s~active) + FactorX(s~A) k_xa_act + + ## AMPLIFICATION Phase (Thrombin Burst) + # Initial Xa produces some Thrombin + FactorX(s~A) + Prothrombin(s~U) -> FactorX(s~A) + Thrombin(s~A) 0.5 + + # Feedback: Thrombin activates Factor V and VIII (Amplification) + Thrombin(s~A) + FactorV(s~U) -> Thrombin(s~A) + FactorV(s~A) k_v_act + FactorX(s~A) + FactorV(s~U) -> FactorX(s~A) + FactorV(s~A) k_v_act # Added starter_act + + # Va and Xa form Prothrombinase complex + FactorV(s~A,b) + FactorX(s~A,b) <-> FactorV(s~A,b!1).FactorX(s~A,b!1) k_complex_va_xa,0.5 + + # Burst: Prothrombinase produces Thrombin rapidly + FactorV(s~A,b!1).FactorX(s~A,b!1) + Prothrombin(s~U) -> FactorV(s~A,b!1).FactorX(s~A,b!1) + Thrombin(s~A) k_pt_burst + + ## PROPAGATION Phase + # Thrombin converts Fibrinogen to Fibrin (TotalRate: aggregate-dependent) + # Wildcard: Enabled if thrombin is active (!?) + Fibrinogen(s~S) -> Fibrinogen(s~F) k_fibrin_synth * (Active_Thrombin/(Km_clot + Active_Thrombin)) TotalRate + + ## TERMINATION Phase (Inhibition) + # Antithrombin neutralizing active Thrombin + # Replaced DeleteMolecules with complex formation to track Inhibitor_Work + Thrombin(s~A) + AT(b) -> Thrombin(s~A!1).AT(b!1) k_at_inh + Thrombin(s~A!1).AT(b!1) -> 0 k_deg_complex + + # Factor clearing + FactorX(s~A) -> FactorX(s~U) 0.01 + FactorV(s~A) -> FactorV(s~U) 0.01 +end reaction rules + +begin actions + generate_network({overwrite=>1}) + simulate({method=>"ode",t_end=>20,n_steps=>200}) +end actions +end model diff --git a/tests/v1/bngl_corpus/bmp-signaling.bngl b/tests/v1/bngl_corpus/bmp-signaling.bngl new file mode 100644 index 00000000..c177a400 --- /dev/null +++ b/tests/v1/bngl_corpus/bmp-signaling.bngl @@ -0,0 +1,105 @@ +begin model +begin parameters + # BMP-Smad signaling: Developmental gradient relay. + # Advanced features: exclude_reactants for Smad6 inhibitor logic. + + # Extracellular Gating + k_noggin_bind 1e-3 # Noggin sequesters BMP (Sink) + k_bmp_bind 1.0 # BMP binds Type II R + + # Receptor Assembly + k_complex 2.0 # Assembly of Type II and Type I + k_smad_phos 3.0 # Phosphorylation of R-Smad (Smad1/5/8) + + # Feedback & Regulation + k_smad6_inh 5.0 # Smad6 inhibits Smad1 recruitment + k_smad6_synth 0.5 # Nuclear signal induces Smad6 + + # Transport + k_import 1.0 # pSmad-Smad4 complex entry + k_export 0.2 + + # Initials + BMP_tot 100 + Noggin_init 50 + R1_tot 150 + R2_tot 150 + Smad1_tot 500 + Smad4_tot 300 + Smad6_init 10 +end parameters + +begin molecule types + BMP(r1,r2,b) + Noggin(b) + Receptor1(l,s~U~P) + Receptor2(l,r,s~U~A) + Smad1(r,s~U~P,loc~cyt~nuc) + Smad4(b,loc~cyt~nuc) + Smad6(b) +end molecule types + +begin seed species + BMP(r1,r2,b) BMP_tot + Noggin(b) Noggin_init + Receptor1(l,s~U) R1_tot + Receptor2(l,r,s~U) R2_tot + Smad1(r,s~U,loc~cyt) Smad1_tot + Smad4(b,loc~cyt) Smad4_tot + Smad6(b) Smad6_init +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules SMAD_Output Smad1(loc~nuc) # Transcription pool + Molecules Active_Nuclear_Smad1 Smad1(loc~nuc,s~P) # Feedback drive status + Molecules Surface_Engage Receptor2(s~A) # Active complex metrics + Molecules Noggin_Sink BMP(b!+) # Sequestration status + Molecules Feedback_Brake Smad6() # Inhibitor levels + Molecules Signal_Complex Smad1(s~P!1).Smad4(b!1) # Functional relay +end observables + +begin functions + # Saturable Smad6 synthesis based on nuclear signal + v_feedback() = k_smad6_synth * (Active_Nuclear_Smad1 / (100 + Active_Nuclear_Smad1)) +end functions + +begin reaction rules + ## EXTRACELLULAR GATING + # Noggin sequesters BMP ligand + BMP(r1,r2,b) + Noggin(b) <-> BMP(r1,r2,b!1).Noggin(b!1) k_noggin_bind,0.05 + + ## RECEPTION + # BMP binds Type II R,then recruits Type I + BMP(r2,b) + Receptor2(l,s~U) <-> BMP(r2!1,b).Receptor2(l!1,s~A) k_bmp_bind,0.1 + Receptor2(l!+,r,s~A) + Receptor1(l,s~U) <-> Receptor2(l!+,r!1,s~A).Receptor1(l!1,s~P) k_complex,0.1 + + ## SIGNAL RELAY + # Active complex phosphorylates Smad1 + # exclude_reactants: Smad1 cannot be phosphorylated if Smad6 is present (Simplified) + Receptor1(s~P) + Smad1(r,s~U) -> Receptor1(s~P) + Smad1(r,s~P) k_smad_phos + + # Inhibitor Smad6 interferes (Sequestration of Smad1) + Smad6(b) + Smad1(r,s~U) <-> Smad6(b!1).Smad1(r!1,s~U) k_smad6_inh,0.2 + + ## TRANSPORT + # pSmad1 dimerizes and translocates + Smad1(s~P,loc~cyt) + Smad4(b,loc~cyt) <-> Smad1(s~P!1,loc~cyt).Smad4(b!1,loc~cyt) 2.0,0.1 + Smad1(loc~cyt,s~U)!+ -> Smad1(loc~nuc,s~U)!+ k_import + + ## FEEDBACK INDUCTION + # Nuclear factor induces Smad6 + Smad1(loc~nuc,s~P) -> Smad1(loc~nuc,s~P) + Smad6(b) v_feedback() + + ## RESET + Smad1(r,loc~nuc,s~P) -> Smad1(r,loc~cyt,s~U) 0.1 + Smad1(loc~nuc) -> Smad1(loc~cyt) k_export + Smad6() -> 0 0.05 + Receptor1(s~P) -> Receptor1(s~U) 0.05 +end reaction rules + +begin actions + generate_network({overwrite=>1}) + simulate({method=>"ode",t_end=>200,n_steps=>400}) +end actions +end model diff --git a/tests/v1/bngl_corpus/brusselator-oscillator.bngl b/tests/v1/bngl_corpus/brusselator-oscillator.bngl new file mode 100644 index 00000000..bb9fb4ba --- /dev/null +++ b/tests/v1/bngl_corpus/brusselator-oscillator.bngl @@ -0,0 +1,94 @@ +begin model +begin parameters + # The Brusselator: Auto-catalytic chemical oscillator. + # Advanced features: Saturable kinetics and if function for phase control. + + # Supply and Decay (Saturable/Michaelis-Menten) + v_supply_max 2.0 # Max input of A + Km_a 0.5 # Reservoir threshold + k_decay_max 1.5 # Max output of X + Km_x 1.2 # Exit threshold + + # Oscillation Drive + k_cross 3.0 # B + X -> Y + k_auto 1.0 # 2X + Y -> 3X (The Positive Feedback) + + # Environmental Control + k_temp_mod 1.0 # Base temp effect + v_temp_current 2.5 + A_pool 100 # Reservoir + B_drive 2.5 # Force + X_init 1.0 + Y_init 1.5 +end parameters + +begin molecule types + A() # Feedstock + B() # Driver + X() # Activator + Y() # Inhibitor/Substrate +end molecule types + +begin seed species + A() A_pool + B() B_drive + X() X_init + Y() Y_init +end seed species + +begin observables + # KEY BIOLOGICAL OUTPUTS + Molecules Activator_X X() # Fast variable + Molecules Inhibitor_Y Y() # Slow variable + Molecules Feed_Status A() # Resource levels + Molecules Drive_Force B() # Control parameter + Molecules Ratio_XY X(),Y() # Phase space metrics + Molecules A_Level A() # Function substrate status + Molecules X_Level X() # Function substrate status +end observables + +begin functions + # Saturable supply from reservoir A + v_supply() = v_supply_max * A_Level / (Km_a + A_Level) + + # Saturable decay of activator X + v_decay() = k_decay_max * X_Level / (Km_x + X_Level) + + + v_auto() = k_auto * v_temp_current +end functions + +begin reaction rules + ## ACTIVATION (Input) + # A -> X (Input flux) + 0 -> X() v_supply() + + ## OSCILLATION Core + # B + X -> Y (Conversion to inhibitor) + B() + X() -> Y() k_cross + + # 2X + Y -> 3X (Auto-catalysis with "Temp" modulation) + # ternary rule - standard Brusselator + X() + X() + Y() -> X() + X() + X() v_auto() + + ## DECAY (Output) + # X -> 0 (Clearance) + X() -> 0 v_decay() + + # Slow resource depletion + A() -> 0 0.05 +end reaction rules + +begin actions + generate_network({overwrite=>1}) + # Phase 1: Cool (0-20) + setParameter("v_temp_current", 2.5) + simulate({method=>"ode",t_end=>20,n_steps=>200}) + # Phase 2: Heat Pulse (20-40) + setParameter("v_temp_current", 5.0) + simulate({method=>"ode",t_end=>40,n_steps=>200,continue=>1}) + # Phase 3: Recovery (40-60) + setParameter("v_temp_current", 2.5) + simulate({method=>"ode",t_end=>60,n_steps=>200,continue=>1}) +end actions +end model diff --git a/tests/v1/bngl_corpus/catalysis.bngl b/tests/v1/bngl_corpus/catalysis.bngl new file mode 100644 index 00000000..da1f2701 --- /dev/null +++ b/tests/v1/bngl_corpus/catalysis.bngl @@ -0,0 +1,81 @@ +# Catalysis in energy BNG +# justin.s.hogg@gmail.com, 9 Apr 2013 + +# requires BioNetGen version >= 2.2.4 +version("2.2.4") +# Quantities have units in moles, so set this to Avogadro's Number +setOption("NumberPerQuantityUnit",6.0221e23) + +begin model +begin parameters + # fundamental constants + RT 2.577 # kJ/mol + NA 6.022e23 # /mol + # simulation volume, L + volC 1e-12 + # initial concentrations, mol/L + conc_S_0 1e-6 + conc_kinase_0 10e-9 + conc_pptase_0 10e-9 + conc_ATP_0 1.0e-3 + conc_ADP_0 0.1e-3 + # standard free energy of formation, kJ/mol + Gf_Sp 51.1 + Gf_S_kinase -41.5 + Gf_S_pptase -41.5 + Gf_ATP 51.1 + # baseline activation energy, kJ/mol + Ea0_S_kinase -7.7 + Ea0_S_pptase -7.7 + Ea0_cat_kinase -11.9 + Ea0_cat_pptase 11.9 + # rate distribution parameter, no units + phi 0.5 +end parameters +begin compartments + # generic compartment + C 3 volC +end compartments +begin molecule types + S(e,y~0~P) # substrate with enzyme binding domain and site of phosphorylation + kinase(s) # kinase enzyme + pptase(s) # phosphotase enzyme + ATP() + ADP() +end molecule types +begin species + S(e,y~0)@C conc_S_0*NA*volC + kinase(s)@C conc_kinase_0*NA*volC + pptase(s)@C conc_pptase_0*NA*volC + $ATP()@C conc_ATP_0*NA*volC # ATP concentration held constant + $ADP()@C conc_ADP_0*NA*volC # ADP concentration held constant +end species +begin reaction rules + # binding rules + S(e) + kinase(s) <-> S(e!1).kinase(s!1) Arrhenius(phi,Ea0_S_kinase) + S(e) + pptase(s) <-> S(e!1).pptase(s!1) Arrhenius(phi,Ea0_S_pptase) + # catalysis + S(e!1,y~0).kinase(s!1) + ATP <-> S(e!1,y~P).kinase(s!1) + ADP Arrhenius(phi,Ea0_cat_kinase) + S(e!1,y~P).pptase(s!1) <-> S(e!1,y~0).pptase(s!1) Arrhenius(phi,Ea0_cat_pptase) +end reaction rules +begin energy patterns + S(y~P) Gf_Sp/RT # phosphorylated subtrate + S(e!0).kinase(s!0) Gf_S_kinase/RT # substrate-kinase binding + S(e!0).pptase(s!0) Gf_S_pptase/RT # substrate-pptase binding + ATP() Gf_ATP/RT # ATP energy (relative to ADP) +end energy patterns +begin observables + Molecules Sp S(y~P) + Molecules S_kinase S(e!1).kinase(s!1) + Molecules S_pptase S(e!1).pptase(s!1) + Molecules Stot S() + Molecules kinaseTot kinase() + Molecules pptaseTot pptase() +end observables +end model + +# generate reaction network.. +generate_network({overwrite=>1}) + +# simulate ODE system to steady state.. +simulate({method=>"ode",t_start=>0,t_end=>3600,n_steps=>120,atol=>1e-3,rtol=>1e-7}) diff --git a/tests/v1/bngl_corpus/egg.bngl b/tests/v1/bngl_corpus/egg.bngl new file mode 100644 index 00000000..6128e684 --- /dev/null +++ b/tests/v1/bngl_corpus/egg.bngl @@ -0,0 +1,64 @@ +# a0__FREE changed to 9.99318747e+01 +# a1__FREE changed to 1.00606018e+00 +# a2__FREE changed to -7.52962956e-01 +# b1__FREE changed to -3.08973913e+01 +# b2__FREE changed to 1.26208508e+00 +# c0__FREE changed to 1.39738978e+02 +# c1__FREE changed to -3.91791542e+01 +# c2__FREE changed to -1.56811184e+00 +# d1__FREE changed to -1.21914775e+00 +# d2__FREE changed to 1.28080122e+00 +# End of permute change log +begin model + begin parameters +a0__FREE 9.99318747e+01 +a1__FREE 1.00606018e+00 +a2__FREE -7.52962956e-01 +b1__FREE -3.08973913e+01 +b2__FREE 1.26208508e+00 +c0__FREE 1.39738978e+02 +c1__FREE -3.91791542e+01 +c2__FREE -1.56811184e+00 +d1__FREE -1.21914775e+00 +d2__FREE 1.28080122e+00 + a0 a0__FREE + a1 a1__FREE + a2 a2__FREE + b1 b1__FREE + b2 b2__FREE + c0 c0__FREE + c1 c1__FREE + c2 c2__FREE + d1 d1__FREE + d2 d2__FREE + pi=2*asin(1) + period 180 + m=2*pi/period + end parameters + begin molecule types + t + end molecule types + begin seed species + t 0 + end seed species + begin observables + Species t t + end observables + begin functions + X()=a0\ + +a1*cos(m*1*t)+b1*sin(m*1*t)\ + +a2*cos(m*2*t)+b2*sin(m*2*t) + Y()=c0\ + +c1*cos(m*1*t)+d1*sin(m*1*t)\ + +c2*cos(m*2*t)+d2*sin(m*2*t) + end functions + begin reaction rules + 0->t 1 + end reaction rules +end model +begin actions + generate_network({overwrite=>1}) + simulate({suffix=>"egg",method=>"ode",\ + t_start=>0,t_end=>180,n_steps=>180,\ + print_functions=>1}) +end actions diff --git a/tests/v1/bngl_corpus/elephant_EFA.bngl b/tests/v1/bngl_corpus/elephant_EFA.bngl new file mode 100644 index 00000000..75ca21ff --- /dev/null +++ b/tests/v1/bngl_corpus/elephant_EFA.bngl @@ -0,0 +1,153 @@ +begin model + begin parameters + a0 48.6399 + a1 30.055 + a2 12.748 + a3 2.2283 + a4 -0.85348 + a5 -1.6215 + a6 -0.17618 + a7 -0.048662 + a8 -1.0801 + a9 0.2382 + a10 -0.13038 + a11 0.12645 + a12 -0.22432 + a13 0.88207 + a14 0.22573 + a15 0.53195 + a16 0.32347 + a17 0.63596 + a18 0.16089 + a19 0.069375 + a20 -0.28683 + b0 0 + b1 21.977 + b2 -3.3539 + b3 -3.5048 + b4 -0.50427 + b5 0.21991 + b6 -0.29889 + b7 -0.50247 + b8 -0.82477 + b9 -0.17609 + b10 0.91897 + b11 -0.39103 + b12 0.60965 + b13 -0.9696 + b14 0.98694 + b15 -0.6184 + b16 0.22323 + b17 0.23272 + b18 0.44101 + b19 0.34473 + b20 -0.43427 + c0 36.6228 + c1 16.478 + c2 -0.69587 + c3 -4.4809 + c4 -5.6869 + c5 0.36898 + c6 4.0373 + c7 -0.013559 + c8 0.010601 + c9 0.36869 + c10 -0.89048 + c11 0.07215 + c12 -0.87898 + c13 -0.51409 + c14 0.28394 + c15 0.13376 + c16 0.2394 + c17 -0.0027284 + c18 -0.058979 + c19 0.1159 + c20 -0.43577 + d0 0 + d1 -15.286 + d2 -10.371 + d3 -3.3776 + d4 -1.3665 + d5 -1.7818 + d6 -8.758 + d7 7.0236 + d8 0.44179 + d9 0.41343 + d10 1.5842 + d11 0.56122 + d12 -0.38208 + d13 0.69121 + d14 0.65805 + d15 0.26396 + d16 0.31597 + d17 -0.20032 + d18 -0.4028 + d19 -0.29874 + d20 0.32228 + pi=2*asin(1) + period 464 + m=2*pi/period + end parameters + begin molecule types + t + end molecule types + begin seed species + t 0 + end seed species + begin observables + Species t t + end observables + begin functions + X()=a0\ + +a1*cos(m*1*t)+b1*sin(m*1*t)\ + +a2*cos(m*2*t)+b2*sin(m*2*t)\ + +a3*cos(m*3*t)+b3*sin(m*3*t)\ + +a4*cos(m*4*t)+b4*sin(m*4*t)\ + +a5*cos(m*5*t)+b5*sin(m*5*t)\ + +a6*cos(m*6*t)+b6*sin(m*6*t)\ + +a7*cos(m*7*t)+b7*sin(m*7*t)\ + +a8*cos(m*8*t)+b8*sin(m*8*t)\ + +a9*cos(m*9*t)+b9*sin(m*9*t)\ + +a10*cos(m*10*t)+b10*sin(m*10*t)\ + +a11*cos(m*11*t)+b11*sin(m*11*t)\ + +a12*cos(m*12*t)+b12*sin(m*12*t)\ + +a13*cos(m*13*t)+b13*sin(m*13*t)\ + +a14*cos(m*14*t)+b14*sin(m*14*t)\ + +a15*cos(m*15*t)+b15*sin(m*15*t)\ + +a16*cos(m*16*t)+b16*sin(m*16*t)\ + +a17*cos(m*17*t)+b17*sin(m*17*t)\ + +a18*cos(m*18*t)+b18*sin(m*18*t)\ + +a19*cos(m*19*t)+b19*sin(m*19*t)\ + +a20*cos(m*20*t)+b20*sin(m*20*t) + Y()=c0\ + +c1*cos(m*1*t)+d1*sin(m*1*t)\ + +c2*cos(m*2*t)+d2*sin(m*2*t)\ + +c3*cos(m*3*t)+d3*sin(m*3*t)\ + +c4*cos(m*4*t)+d4*sin(m*4*t)\ + +c5*cos(m*5*t)+d5*sin(m*5*t)\ + +c6*cos(m*6*t)+d6*sin(m*6*t)\ + +c7*cos(m*7*t)+d7*sin(m*7*t)\ + +c8*cos(m*8*t)+d8*sin(m*8*t)\ + +c9*cos(m*9*t)+d9*sin(m*9*t)\ + +c10*cos(m*10*t)+d10*sin(m*10*t)\ + +c11*cos(m*11*t)+d11*sin(m*11*t)\ + +c12*cos(m*12*t)+d12*sin(m*12*t)\ + +c13*cos(m*13*t)+d13*sin(m*13*t)\ + +c14*cos(m*14*t)+d14*sin(m*14*t)\ + +c15*cos(m*15*t)+d15*sin(m*15*t)\ + +c16*cos(m*16*t)+d16*sin(m*16*t)\ + +c17*cos(m*17*t)+d17*sin(m*17*t)\ + +c18*cos(m*18*t)+d18*sin(m*18*t)\ + +c19*cos(m*19*t)+d19*sin(m*19*t)\ + +c20*cos(m*20*t)+d20*sin(m*20*t) + end functions + begin reaction rules + 0->t 1 + end reaction rules +end model +begin actions + generate_network({overwrite=>1}) + simulate({suffix=>"elephant",method=>"ode",\ + t_start=>0,t_end=>464,n_steps=>464,\ + print_functions=>1}) +end actions diff --git a/tests/v1/bngl_corpus/energy_transport_pump.bngl b/tests/v1/bngl_corpus/energy_transport_pump.bngl new file mode 100644 index 00000000..8d723821 --- /dev/null +++ b/tests/v1/bngl_corpus/energy_transport_pump.bngl @@ -0,0 +1,85 @@ +# Model: energy_transport_pump.bngl +# Description: Models an active transport pump (T) that moves substrate A against its chemical potential gradient. +# The process is driven by the hydrolysis of ATP. +# Overall reaction: A(out) + ATP -> A(in) + ADP + Pi +# Illustrates how coupling to a high-energy reaction can drive an unfavorable transition. + +begin model + +begin parameters + # --- Transport & Chemical Energies --- + # A(out) is the ground/low energy state + G_A_out 0 + # A(in) is the high energy state (concentration gradient or electrical potential) + G_A_in 5.0 + + # ATP Hydrolysis Energies (High Energy) + G_ATP 10.0 + G_ADP 0.0 + G_Pi 0.0 + + # Activation barrier for the coupled transport cycle + Ea_transport 2.0 + + # System + RT 1.0 + phi 0.5 + + # --- Concentrations --- + Atot 100 + ATP_conc 500 # High ATP pool + ADP_conc 10 + Pi_conc 10 +end parameters + +begin molecule types + A(loc~in~out) # Substrate A with location state + ATP() + ADP() + Pi() + T(a) # Transporter Enzyme +end molecule types + +begin seed species + A(loc~out) Atot # Start with all A outside + ATP() ATP_conc + ADP() ADP_conc + Pi() Pi_conc + T(a) 10 # Limited number of transporters +end seed species + +begin observables + Molecules A_in A(loc~in) + Molecules A_out A(loc~out) + Molecules Energy_Source ATP() +end observables + +begin energy patterns + # Define energies for reactants and products + A(loc~out) G_A_out + A(loc~in) G_A_in + ATP() G_ATP + ADP() G_ADP + Pi() G_Pi + # Transporter is catalyst, G_T cancels out on both sides +end energy patterns + +begin reaction rules + # Coupled Transport Rule + # Reactants: T + A(out) + ATP + # Products: T + A(in) + ADP + Pi + # Delta G = (G_A_in + G_ADP + G_Pi) - (G_A_out + G_ATP) + # = (5 + 0 + 0) - (0 + 10) = -5.0 + # Since Delta G is negative, the forward reaction is spontaneous. + T(a) + A(loc~out) + ATP() <-> T(a) + A(loc~in) + ADP() + Pi() Arrhenius(phi, Ea_transport) + + # Passive Leak Channel (Backflow) + # A(in) -> A(out) is favorable (Delta G = -5), but we put a high barrier to minimize it. + A(loc~in) <-> A(loc~out) Arrhenius(phi, 8.0) +end reaction rules + +end model + +## Actions ## +generate_network({overwrite=>1}) +simulate({method=>"ode", t_end=>20, n_steps=>200}) diff --git a/tests/v1/bngl_corpus/example1.bngl b/tests/v1/bngl_corpus/example1.bngl new file mode 100644 index 00000000..ac8a1c7c --- /dev/null +++ b/tests/v1/bngl_corpus/example1.bngl @@ -0,0 +1,58 @@ +# Example file for BNG2 tutorial. +# All text following the occurence off '#' character in a line is ignored. +# Written by James R. Faeder +# Theoretical Biology and Biophysics Group +# Los Alamos National Laboratory +# faeder@lanl.gov +# 10/28/2005 +# revised 6/9/2006 + +version("2.0.23"); + +begin parameters + 1 L0 1 + 2 R0 1 + 3 A0 5 + 4 kp1 0.5 + 5 km1 0.1 + 6 kp2 1.1 + 7 km2 0.1 + 8 p1 10 + 9 d1 5 + 10 kpA 1e1 + 11 kmA 0.02 +end parameters + +begin species + 1 L(r) L0 + 2 R(l,d,Y~U) R0 + 3 A(SH2) A0 +end species + +begin observables + Molecules R_dim R(d!+) + Molecules R_phos R(Y~P!?) + Molecules A_R A(SH2!1).R(Y~P!1) +end observables + +begin reaction rules + 1 L(r) + R(l,d) <-> L(r!1).R(l!1,d) kp1, km1 + 2 R(l!+,d) + R(l!+,d) <-> R(l!+,d!2).R(l!+,d!2) kp2, km2 + 3 R(d!+,Y~U) -> R(d!+,Y~P) p1 + 4 R(Y~P) -> R(Y~U) d1 + 5 R(Y~P) + A(SH2) <-> R(Y~P!1).A(SH2!1) kpA, kmA +end reaction rules + +# Call with no arguments +generate_network(); + +# Call with a single parameter +#generate_network({max_iter=>2}); + +# Call with a hash valued parameter +#generate_network({max_stoich=>{A=>1}}); + +simulate_ode({t_end=>50,n_steps=>20}); + +# Print concentratons at unevenly spaced times (array-valued parameter) +#simulate_ode({sample_times=>[1,10,100]}); diff --git a/tests/v1/bngl_corpus/genetic_bistability_energy.bngl b/tests/v1/bngl_corpus/genetic_bistability_energy.bngl new file mode 100644 index 00000000..d34468f1 --- /dev/null +++ b/tests/v1/bngl_corpus/genetic_bistability_energy.bngl @@ -0,0 +1,91 @@ +# Model: genetic_bistability_energy.bngl +# Description: Models a bistable genetic toggle switch where two genes (A and B) repress each other. +# Uses Energy Patterns to model the cooperative repression (high affinity binding). +# Demonstrates how to build stability landscapes. + +begin model + +begin parameters + # Energy + G_Gene 0 + G_Protein 0 + + # Binding Energy (Repressor to Operator) + # Strong binding + G_bind -5.0 + G_coop -2.0 # Dimerization/Cooperativity + + RT 1.0 + phi 0.5 + Ea_bind 1.0 + + k_synth 10.0 + k_deg 0.1 +end parameters + +begin molecule types + GeneA(op) # Operator site + GeneB(op) + ProtA(r,d) # Repressor domain, Dimer domain + ProtB(r,d) +end molecule types + +begin seed species + # Start in unstable middle state + GeneA(op) 1 + GeneB(op) 1 + ProtA(r,d) 10 + ProtB(r,d) 10 +end seed species + +begin observables + Molecules A ProtA() + Molecules B ProtB() + Molecules GeneA_Free GeneA(op) + Molecules GeneB_Free GeneB(op) +end observables + +begin energy patterns + GeneA(op) G_Gene + GeneB(op) G_Gene + ProtA(r,d) G_Protein + ProtB(r,d) G_Protein + + # Specific Interactions + # ProtA binds GeneB + GeneB(op!1).ProtA(r!1) G_bind + # ProtB binds GeneA + GeneA(op!1).ProtB(r!1) G_bind + + # Dimerization (Cooperative) + ProtA(d!1).ProtA(d!1) G_coop + ProtB(d!1).ProtB(d!1) G_coop +end energy patterns + +begin reaction rules + # 1. Synthesis + # Only free genes express + GeneA(op) -> GeneA(op) + ProtA(r,d) k_synth + GeneB(op) -> GeneB(op) + ProtB(r,d) k_synth + + # 2. Binding (Repression) + # Rates governed by Energy Patterns + # ProtA binds GeneB + GeneB(op) + ProtA(r) <-> GeneB(op!1).ProtA(r!1) Arrhenius(phi, Ea_bind) + # ProtB binds GeneA + GeneA(op) + ProtB(r) <-> GeneA(op!1).ProtB(r!1) Arrhenius(phi, Ea_bind) + + # 3. Degradation + ProtA() -> 0 k_deg + ProtB() -> 0 k_deg + + # 4. Dimerization (Optional, adds non-linearity) + ProtA(d) + ProtA(d) <-> ProtA(d!1).ProtA(d!1) Arrhenius(phi, 1.0) + ProtB(d) + ProtB(d) <-> ProtB(d!1).ProtB(d!1) Arrhenius(phi, 1.0) +end reaction rules + +end model + +## Actions ## +generate_network({overwrite=>1}) +simulate({method=>"ode", t_end=>100, n_steps=>200}) diff --git a/tests/v1/bngl_corpus/golden.json b/tests/v1/bngl_corpus/golden.json new file mode 100644 index 00000000..c702e21a --- /dev/null +++ b/tests/v1/bngl_corpus/golden.json @@ -0,0 +1,1429 @@ +{ + "An_2009.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "A20", + "A20mRNA", + "Administer", + "CD14", + "DNA", + "IRAK1", + "IRAK4", + "IRAKM", + "IkB", + "IkBmRNA", + "Ikk_Complex", + "LPS", + "MAL", + "MD2", + "MyD88", + "MyD88s", + "NFkB", + "Proteasome26s", + "RP1", + "SARM", + "TAK1", + "TLR4", + "TNF", + "TNFmRNA", + "TRAF4", + "TRAF6", + "TRAM", + "TRIF", + "Tollip", + "Trash", + "iNOSmRNA" + ], + "observables": [ + "A20", + "A20_NFkB_DNA", + "Activated_Ikk_complex", + "Activated_TAK1", + "IkB_Degraded", + "IkB_Prot26s", + "IkB_active", + "IkBmRNA_Off", + "Inactive_Cyto_NFkB", + "NFkB_Active_Cyto", + "NFkB_Active_Nucleus", + "NFkB_DNA_IkB", + "NFkB_Inactive", + "NonBoundNonPhos_IkB", + "Phos_IkB_NFkB", + "TNF", + "TNF_NFkB_DNA", + "TNFmRNA_Off", + "Unbound_Cyto_NFkB" + ], + "parameters": [ + "A20_Degrade", + "A20_IkkAct_Deactivate", + "A20_Init", + "A20_MyD88IRAK1TRAF6_Degrade", + "A20_Preconditioned", + "A20_TRAF6TRIFRP1_Degrade", + "A20_TRAF6_Bind", + "A20_TRAF6_Unbind", + "A20_Transcription_Execute", + "A20_Translation_Execute", + "CD14_Init", + "CD14_MD2_Bind", + "CD14_MD2_Unbind", + "CD14_TLR4_Bind", + "CD14_TLR4_Unbind", + "DNA", + "IRAK1_IRAK4_Bind", + "IRAK1_IRAK4_Unbind", + "IRAK1_Init", + "IRAK4_Init", + "IkB_DegradeNFkB", + "IkB_Init", + "IkB_Proteasome23_Degrade", + "IkB_Transcription_Execute", + "IkB_Translation_Execute", + "Ikk_Complex_Init", + "Ikk_Deactivation", + "Ikk_Degradation_Rate", + "Ikk_complex_IkB_Phos", + "LPS_CD14_Bind", + "LPS_CD14_Unbind", + "LPS_Init", + "LPS_MD2_Bind", + "LPS_MD2_Unbind", + "LPS_TLR4_Bind", + "LPS_TLR4_Unbind", + "MAL_Init", + "MD2_Init", + "MD2_TLR4_Bind", + "MD2_TLR4_Unbind", + "MyD88IRAK1TRAF6_TAK1_Activate", + "MyD88IRAK1_TRAF6_B_Unbind", + "MyD88IRAK1_TRAF6_Bind", + "MyD88_IRAK1_Bind", + "MyD88_IRAK1_Unbind", + "MyD88_IRAK4_Bind", + "MyD88_IRAK4_Unbind", + "MyD88_Init", + "NFkB_DNA_A20_Bind", + "NFkB_DNA_A20_Unbind", + "NFkB_DNA_IkB_Bind", + "NFkB_DNA_IkB_Unbind", + "NFkB_DNA_TNF_Bind", + "NFkB_DNA_TNF_Unbind", + "NFkB_Degredation", + "NFkB_IkB_Bind", + "NFkB_IkB_Unbind", + "NFkB_Inactive_Cytoplasm", + "NFkB_Translocation_Nucleus", + "Proteasome23_Init", + "RP1_Init", + "RP1_TRAF6_Bind", + "RP1_TRAF6_Unbind", + "RP1_TRIF_Bind", + "RP1_TRIF_Unbind", + "TAK1_Deactivation", + "TAK1_Degradation", + "TAK1_Ikk_Complex_Activate", + "TAK1_Init", + "TLR4MAL_MyD88_Bind", + "TLR4MAL_MyD88_Unbind", + "TLR4TRAM_TRIF_Bind", + "TLR4TRAM_TRIF_Unbind", + "TLR4_Complex_Dimer_Bind", + "TLR4_Complex_Dimer_Unbind", + "TLR4_MAL_Bind", + "TLR4_MAL_Unbind", + "TLR4_TRAM_Bind", + "TLR4_TRAM_Unbind", + "TLR_Init", + "TNF_Degrade", + "TNF_Transcription_Execute", + "TNF_Translation_Execute", + "TRAF6TRIF_TAK1_Activate", + "TRAF6_Init", + "TRAF6_MyD88IRAK1_Bind", + "TRAF6_MyD88IRAK1_Unbind", + "TRAF6_TRIF_Bind", + "TRAF6_TRIF_Unbind", + "TRAM_Init", + "TRIF_Init", + "TRIF_TRAF6_Bind", + "TRIF_TRAF6_Unbind", + "p50_Init", + "p65_Init", + "p65_p50_Bind", + "p65_p50_Unbind" + ], + "seed_composition": [ + [ + "A20" + ], + [ + "CD14" + ], + [ + "DNA" + ], + [ + "IRAK1" + ], + [ + "IRAK4" + ], + [ + "IkB" + ], + [ + "Ikk_Complex" + ], + [ + "LPS" + ], + [ + "MAL" + ], + [ + "MD2" + ], + [ + "MyD88" + ], + [ + "NFkB" + ], + [ + "Proteasome26s" + ], + [ + "RP1" + ], + [ + "TAK1" + ], + [ + "TLR4" + ], + [ + "TRAF6" + ], + [ + "TRAM" + ], + [ + "TRIF" + ] + ] + }, + "Barua_2009__PATCHED.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "J", + "S" + ], + "observables": [ + "JS", + "JSS", + "JSSJ", + "J_active", + "J_inactive", + "J_mono" + ], + "parameters": [ + "Jtot", + "Stot", + "koff_SH2", + "koff_dimer", + "kon_SH2", + "kon_dimer", + "kphos_fast", + "kphos_slow" + ], + "seed_composition": [ + [ + "J" + ], + [ + "S" + ] + ] + }, + "Chattaraj_2021.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "NWASP", + "Nck", + "Nephrin" + ], + "observables": [ + "cluster_nck_nw", + "cluster_neph_nck_nw", + "free_NWASP", + "free_Nck", + "free_Nephrin", + "fully_bound_NWASP", + "fully_bound_Nck", + "fully_bound_Nephrin", + "tot_NWASP", + "tot_Nck", + "tot_Nephrin" + ], + "parameters": [ + "kd_12", + "kd_23", + "koff_12", + "koff_23", + "kon_12", + "kon_23" + ], + "seed_composition": [ + [ + "NWASP" + ], + [ + "Nck" + ], + [ + "Nephrin" + ] + ] + }, + "LR.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "L", + "R" + ], + "observables": [ + "Bound", + "FreeR" + ], + "parameters": [ + "L0", + "NaV", + "R0", + "Vcell", + "Vec", + "km1", + "kp1", + "lig_conc" + ], + "seed_composition": [ + [ + "L" + ], + [ + "R" + ] + ] + }, + "LRR.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "L", + "R" + ], + "observables": [ + "Bound", + "Dimers", + "FreeR" + ], + "parameters": [ + "Acell", + "L0", + "NaV", + "R0", + "Vcell", + "Vec", + "Vpm", + "d_pm", + "km1", + "kp1", + "lig_conc" + ], + "seed_composition": [ + [ + "L" + ], + [ + "R" + ] + ] + }, + "Motivating_example_cBNGL.bngl": { + "compartments": [ + "CP", + "EC", + "EM", + "EN", + "NM", + "NU", + "PM" + ], + "functions": [], + "molecule_types": [ + "DNA", + "Im", + "L", + "NP", + "P1", + "P2", + "R", + "Sink", + "TF", + "mRNA1", + "mRNA2" + ], + "observables": [ + "Bound_prom1", + "Bound_prom2", + "Catalytic_R", + "Catalytic_TF", + "CountSink", + "Im_CP", + "Im_Cargo_NP", + "Im_NU", + "L_Bound_EM", + "L_Bound_PM", + "L_Dimers_EC", + "L_Dimers_EM", + "L_Dimers_EN", + "L_Dimers_PM", + "P1_CP", + "P1_NU", + "P1_NU_dna", + "P1_NU_free", + "Phos_TF", + "R_Dimers_EM", + "R_Dimers_PM", + "TF_Dimer_CP", + "TF_Dimer_NU", + "Tot_DNA", + "Tot_Im", + "Tot_L", + "Tot_NP", + "Tot_P1", + "Tot_P2", + "Tot_R", + "Tot_TF", + "Tot_mRNA1", + "Tot_mRNA2" + ], + "parameters": [ + "DNA0", + "Im0", + "L0", + "NP0", + "R0", + "TF0", + "eff_width", + "k_Im_bind_CP", + "k_Im_bind_NU", + "k_Im_cross_NP", + "k_Im_enters_NP", + "k_Im_exits_NP", + "k_Im_unbind_CP", + "k_Im_unbind_NU", + "k_P_deg", + "k_R_dephos", + "k_R_endo", + "k_R_transphos", + "k_TF_dephos", + "k_TF_transphos", + "k_mRNA_deg", + "k_mRNA_to_CP", + "k_recycle", + "k_transcribe", + "k_translate", + "km_LL", + "km_LR", + "km_P1_p2", + "km_R_TF", + "km_R_TFp", + "km_TF_TF", + "km_TF_p1", + "kp_LL", + "kp_LR", + "kp_P1_p2", + "kp_R_TF", + "kp_R_TFp", + "kp_TF_TF", + "kp_TF_p1", + "nEndo", + "sa_EM", + "sa_NM", + "sa_PM", + "vol_CP", + "vol_EC", + "vol_EN", + "vol_NU" + ], + "seed_composition": [ + [ + "DNA" + ], + [ + "Im" + ], + [ + "L" + ], + [ + "NP" + ], + [ + "R" + ], + [ + "Sink" + ], + [ + "TF" + ] + ] + }, + "Ras_bistability_v2.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "CRAF", + "GDP", + "GTP", + "KRAS", + "KRASG12D", + "RASA1", + "RTK", + "SOS1" + ], + "observables": [ + "ActiveEffector", + "CRAFtot", + "Gtot", + "KRAStot", + "RASA1tot", + "RTKbndGAP", + "RTKbndGEF", + "RTKtot", + "RasGDP", + "RasGTP", + "SOS1tot" + ], + "parameters": [ + "CRAF_num", + "GDP_num", + "GTP_num", + "K1", + "K5", + "KRAS_WT_num", + "KRAS_mut_num", + "KRAS_tot_num", + "NA", + "RASA1_num", + "RTK_num", + "SOS1_num", + "Vcyto", + "Vpm", + "f", + "fmut", + "kaEff", + "kaGAP_cyto", + "kaGAP_pm", + "kaGEF_cyto", + "kaGEF_pm", + "kcat", + "kdEff", + "kdGAP", + "kdGEF", + "kf_cyto", + "kf_pm", + "khyd", + "km1", + "km2", + "km3", + "km4", + "km5", + "km6", + "km7", + "kp1", + "kp2_cyto", + "kp2_pm", + "kp3_cyto", + "kp3_pm", + "kp4", + "kp5", + "kp6", + "kp7_cyto", + "kp7_pm", + "kr" + ], + "seed_composition": [ + [ + "CRAF" + ], + [ + "GDP" + ], + [ + "GTP" + ], + [ + "KRAS" + ], + [ + "KRAS" + ], + [ + "KRAS" + ], + [ + "KRASG12D" + ], + [ + "KRASG12D" + ], + [ + "KRASG12D" + ], + [ + "RASA1" + ], + [ + "RTK" + ], + [ + "SOS1" + ] + ] + }, + "Rule_based_egfr_tutorial.bngl": { + "compartments": [ + "c0" + ], + "functions": [], + "molecule_types": [ + "EGF", + "EGFR", + "Grb2", + "Shc" + ], + "observables": [ + "Dimers", + "Dimers_s", + "O0_EGFR_tot", + "O0_EGF_tot", + "O0_Grb2_tot", + "O0_Shc_tot", + "Y1", + "Y2", + "Y_total" + ], + "parameters": [], + "seed_composition": [ + [ + "EGF" + ], + [ + "EGFR" + ], + [ + "Grb2" + ], + [ + "Shc" + ], + [ + "Shc" + ] + ] + }, + "akt-signaling.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "AKT", + "GrowthFactor", + "PI3K", + "RTK", + "S6K", + "mTORC1", + "mTORC2" + ], + "observables": [ + "Active_PI3K", + "Active_RTK", + "Active_mTORC1", + "Double_pAkt", + "pAkt_S473", + "pAkt_T308" + ], + "parameters": [ + "AKT_tot", + "GF_tot", + "PI3K_tot", + "RTK_tot", + "S6K_tot", + "k_akt_s473", + "k_akt_t308", + "k_bind", + "k_fb", + "k_mtorc1", + "k_pi3k", + "k_reset", + "k_rtk_act", + "k_s6k_act", + "k_unbind", + "mTORC1_tot", + "mTORC2_tot" + ], + "seed_composition": [ + [ + "AKT" + ], + [ + "GrowthFactor" + ], + [ + "PI3K" + ], + [ + "RTK" + ], + [ + "S6K" + ], + [ + "mTORC1" + ], + [ + "mTORC2" + ] + ] + }, + "apoptosis-cascade.bngl": { + "compartments": [], + "functions": [ + "v_momp" + ], + "molecule_types": [ + "Apaf1", + "Bid", + "Caspase3", + "Caspase8", + "DeathLigand", + "Mito", + "SMAC", + "XIAP" + ], + "observables": [ + "Active_Bid", + "CytoC_Released", + "Global_Death", + "Inhibitor_Sequest", + "Mitochondria_Gate", + "Mitochondrial_Fail" + ], + "parameters": [ + "Apaf1_tot", + "Bid_tot", + "C3_tot", + "C8_tot", + "Km_tbid", + "MOMP_sites", + "SMAC_tot", + "XIAP_tot", + "k_apaf_act", + "k_bid_trunc", + "k_c3_act", + "k_ligand", + "k_momp_max", + "k_smac_rec", + "k_xiap_bind", + "n_hill" + ], + "seed_composition": [ + [ + "Apaf1" + ], + [ + "Bid" + ], + [ + "Caspase3" + ], + [ + "Caspase8" + ], + [ + "DeathLigand" + ], + [ + "Mito" + ], + [ + "SMAC" + ], + [ + "XIAP" + ] + ] + }, + "bcr-signaling.bngl": { + "compartments": [], + "functions": [ + "v_ca" + ], + "molecule_types": [ + "Antigen", + "BCR", + "CD22", + "Calcium", + "PLCg2", + "SHP1", + "Syk" + ], + "observables": [ + "Active_PLCg2", + "Active_Relay", + "Calcium_Signal", + "Effector_Load", + "Negative_Brake", + "Sync_Antigen" + ], + "parameters": [ + "Antigen_tot", + "BCR_tot", + "CD22_tot", + "Ca_cyt", + "Km_ca", + "PLCg2_tot", + "SHP1_tot", + "Syk_tot", + "k_bind_ag", + "k_ca_max", + "k_dephos_bcr", + "k_endo", + "k_lyn_phos", + "k_plcg_act", + "k_shp_rec", + "k_syk_act", + "n_hill" + ], + "seed_composition": [ + [ + "Antigen" + ], + [ + "BCR" + ], + [ + "CD22" + ], + [ + "Calcium" + ], + [ + "PLCg2" + ], + [ + "SHP1" + ], + [ + "Syk" + ] + ] + }, + "blood-coagulation-thrombin.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "AT", + "FactorV", + "FactorX", + "Fibrinogen", + "Prothrombin", + "TF", + "Thrombin" + ], + "observables": [ + "Active_Thrombin", + "Clot_Mass", + "Feedback_Factors", + "Inhibitor_Work", + "Proth_Complex" + ], + "parameters": [ + "ATIII_tot", + "FactorV_tot", + "FactorX_tot", + "Fibrino_tot", + "Km_clot", + "Prothrom_tot", + "TF_sites", + "k_at_inh", + "k_complex_va_xa", + "k_deg_complex", + "k_fibrin_synth", + "k_pt_burst", + "k_tf_vii", + "k_tpi_inh", + "k_v_act", + "k_viii_act", + "k_xa_act" + ], + "seed_composition": [ + [ + "AT" + ], + [ + "FactorV" + ], + [ + "FactorX" + ], + [ + "Fibrinogen" + ], + [ + "Prothrombin" + ], + [ + "TF" + ] + ] + }, + "bmp-signaling.bngl": { + "compartments": [], + "functions": [ + "v_feedback" + ], + "molecule_types": [ + "BMP", + "Noggin", + "Receptor1", + "Receptor2", + "Smad1", + "Smad4", + "Smad6" + ], + "observables": [ + "Active_Nuclear_Smad1", + "Feedback_Brake", + "Noggin_Sink", + "SMAD_Output", + "Signal_Complex", + "Surface_Engage" + ], + "parameters": [ + "BMP_tot", + "Noggin_init", + "R1_tot", + "R2_tot", + "Smad1_tot", + "Smad4_tot", + "Smad6_init", + "k_bmp_bind", + "k_complex", + "k_export", + "k_import", + "k_noggin_bind", + "k_smad6_inh", + "k_smad6_synth", + "k_smad_phos" + ], + "seed_composition": [ + [ + "BMP" + ], + [ + "Noggin" + ], + [ + "Receptor1" + ], + [ + "Receptor2" + ], + [ + "Smad1" + ], + [ + "Smad4" + ], + [ + "Smad6" + ] + ] + }, + "brusselator-oscillator.bngl": { + "compartments": [], + "functions": [ + "v_auto", + "v_decay", + "v_supply" + ], + "molecule_types": [ + "A", + "B", + "X", + "Y" + ], + "observables": [ + "A_Level", + "Activator_X", + "Drive_Force", + "Feed_Status", + "Inhibitor_Y", + "Ratio_XY", + "X_Level" + ], + "parameters": [ + "A_pool", + "B_drive", + "Km_a", + "Km_x", + "X_init", + "Y_init", + "k_auto", + "k_cross", + "k_decay_max", + "k_temp_mod", + "v_supply_max", + "v_temp_current" + ], + "seed_composition": [ + [ + "A" + ], + [ + "B" + ], + [ + "X" + ], + [ + "Y" + ] + ] + }, + "catalysis.bngl": { + "compartments": [ + "C" + ], + "functions": [], + "molecule_types": [ + "ADP", + "ATP", + "S", + "kinase", + "pptase" + ], + "observables": [ + "S_kinase", + "S_pptase", + "Sp", + "Stot", + "kinaseTot", + "pptaseTot" + ], + "parameters": [ + "Ea0_S_kinase", + "Ea0_S_pptase", + "Ea0_cat_kinase", + "Ea0_cat_pptase", + "Gf_ATP", + "Gf_S_kinase", + "Gf_S_pptase", + "Gf_Sp", + "NA", + "RT", + "conc_ADP_0", + "conc_ATP_0", + "conc_S_0", + "conc_kinase_0", + "conc_pptase_0", + "phi", + "volC" + ], + "seed_composition": [ + [ + "ADP" + ], + [ + "ATP" + ], + [ + "S" + ], + [ + "kinase" + ], + [ + "pptase" + ] + ] + }, + "egg.bngl": { + "compartments": [], + "functions": [ + "X", + "Y" + ], + "molecule_types": [ + "t" + ], + "observables": [ + "t" + ], + "parameters": [ + "a0", + "a0__FREE", + "a1", + "a1__FREE", + "a2", + "a2__FREE", + "b1", + "b1__FREE", + "b2", + "b2__FREE", + "c0", + "c0__FREE", + "c1", + "c1__FREE", + "c2", + "c2__FREE", + "d1", + "d1__FREE", + "d2", + "d2__FREE", + "m", + "period", + "pi" + ], + "seed_composition": [ + [ + "t" + ] + ] + }, + "elephant_EFA.bngl": { + "compartments": [], + "functions": [ + "X", + "Y" + ], + "molecule_types": [ + "t" + ], + "observables": [ + "t" + ], + "parameters": [ + "a0", + "a1", + "a10", + "a11", + "a12", + "a13", + "a14", + "a15", + "a16", + "a17", + "a18", + "a19", + "a2", + "a20", + "a3", + "a4", + "a5", + "a6", + "a7", + "a8", + "a9", + "b0", + "b1", + "b10", + "b11", + "b12", + "b13", + "b14", + "b15", + "b16", + "b17", + "b18", + "b19", + "b2", + "b20", + "b3", + "b4", + "b5", + "b6", + "b7", + "b8", + "b9", + "c0", + "c1", + "c10", + "c11", + "c12", + "c13", + "c14", + "c15", + "c16", + "c17", + "c18", + "c19", + "c2", + "c20", + "c3", + "c4", + "c5", + "c6", + "c7", + "c8", + "c9", + "d0", + "d1", + "d10", + "d11", + "d12", + "d13", + "d14", + "d15", + "d16", + "d17", + "d18", + "d19", + "d2", + "d20", + "d3", + "d4", + "d5", + "d6", + "d7", + "d8", + "d9", + "m", + "period", + "pi" + ], + "seed_composition": [ + [ + "t" + ] + ] + }, + "energy_transport_pump.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "A", + "ADP", + "ATP", + "Pi", + "T" + ], + "observables": [ + "A_in", + "A_out", + "Energy_Source" + ], + "parameters": [ + "ADP_conc", + "ATP_conc", + "Atot", + "Ea_transport", + "G_ADP", + "G_ATP", + "G_A_in", + "G_A_out", + "G_Pi", + "Pi_conc", + "RT", + "phi" + ], + "seed_composition": [ + [ + "A" + ], + [ + "ADP" + ], + [ + "ATP" + ], + [ + "Pi" + ], + [ + "T" + ] + ] + }, + "example1.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [], + "observables": [ + "A_R", + "R_dim", + "R_phos" + ], + "parameters": [ + "A0", + "L0", + "R0", + "d1", + "km1", + "km2", + "kmA", + "kp1", + "kp2", + "kpA", + "p1" + ], + "seed_composition": [ + [ + "A" + ], + [ + "L" + ], + [ + "R" + ] + ] + }, + "genetic_bistability_energy.bngl": { + "compartments": [], + "functions": [], + "molecule_types": [ + "GeneA", + "GeneB", + "ProtA", + "ProtB" + ], + "observables": [ + "A", + "B", + "GeneA_Free", + "GeneB_Free" + ], + "parameters": [ + "Ea_bind", + "G_Gene", + "G_Protein", + "G_bind", + "G_coop", + "RT", + "k_deg", + "k_synth", + "phi" + ], + "seed_composition": [ + [ + "GeneA" + ], + [ + "GeneB" + ], + [ + "ProtA" + ], + [ + "ProtB" + ] + ] + }, + "immob_equiv_lig_sites.bngl": { + "compartments": [], + "functions": [ + "avg_agg_size_R1", + "avg_agg_size_R2", + "large_T2", + "mobility_factor_T2", + "twobigprod_T2a", + "twobigprod_T2b" + ], + "molecule_types": [ + "L1", + "L2", + "R1", + "R2" + ], + "observables": [ + "L1free", + "L2free", + "R1_naggs", + "R1free", + "R1tot", + "R2_naggs", + "R2free", + "R2tot", + "T2aL1", + "T2aL2", + "T2aL3", + "T2aL4", + "T2aL5", + "T2aL6", + "T2aL7", + "T2aL8", + "T2aS1", + "T2aS2", + "T2aS3", + "T2aS4", + "T2bL1", + "T2bL2", + "T2bL3", + "T2bL4", + "T2bL5", + "T2bL6", + "T2bL7", + "T2bL8", + "T2bS1", + "T2bS2", + "T2bS3", + "T2bS4", + "monomericR1", + "monomericR2" + ], + "parameters": [ + "Lconc", + "Lconc_nM", + "LcopyNum", + "NA", + "RcopyNum", + "T2", + "T3", + "Vecf", + "delta", + "epsilon", + "f", + "kf", + "kmx", + "kpx", + "kr", + "rho", + "sigma", + "valL", + "valR" + ], + "seed_composition": [ + [ + "L1" + ], + [ + "L2" + ], + [ + "R1" + ], + [ + "R2" + ] + ] + } +} diff --git a/tests/v1/bngl_corpus/immob_equiv_lig_sites.bngl b/tests/v1/bngl_corpus/immob_equiv_lig_sites.bngl new file mode 100644 index 00000000..bfa2b788 --- /dev/null +++ b/tests/v1/bngl_corpus/immob_equiv_lig_sites.bngl @@ -0,0 +1,298 @@ +# filename: immob_equiv_lig_site.bngl + +begin model + +# See diffusional_slowing_v2.bngl for a version of this model that works. +# In the diffusional_slowing_v2.bngl file, ligand sites are not equivalent. +# Here, ligand sites are equivalent. +# Both models are for bivalent ligand-bivalent receptor interaction with +# immobilization of ligand-induced receptor aggregates when aggregates contain +# more than a threshold number (2 or 3) receptors. + +begin parameters + +# fraction of cell to consider in stochastic simulation +f 0.01 # dimensionless, 01 + +Species L1free L1(r,r) +Species R1free R1(l1,l2) +Molecules R1tot R1() +Species monomericR1 R1==1 +Species R1_naggs R1>1 + +Species L2free L2(r,r) +Species R2free R2(a,b) +Molecules R2tot R2() +Species monomericR2 R2==1 +Species R2_naggs R2>1 + +#Species L3free L3(r,r) +#Species R3free R3(a,b) +#Molecules R3tot R3() +#Species monomericR3 R3==1 +#Species R3_naggs R3>1 + +# observables used in (local) indicator functions for the T2=2 and T3=3 cases + +######## +# T2=2 # +######## + +#Molecules P_T2_aR_R3L2_aaa R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4) +#Molecules P_T2_aR_R3L2_aab R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4) +#Molecules P_T2_aR_R3L2_aba R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4) +#Molecules P_T2_aR_R3L2_abb R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4) + +#Molecules P_T2_La_R3L3_aaa R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_baa R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_aba R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_bba R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_aab R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_bab R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_abb R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +#Molecules P_T2_La_R3L3_bbb R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) + +#Molecules P_T2_Rb_R3L2_aaa R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +#Molecules P_T2_Rb_R3L2_baa R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +#Molecules P_T2_Rb_R3L2_aba R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +#Molecules P_T2_Rb_R3L2_bba R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) + +#Molecules P_T2_bL_R3L3_aaa R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6) +#Molecules P_T2_bL_R3L3_aab R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(b!6) +#Molecules P_T2_bL_R3L3_aba R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6) +#Molecules P_T2_bL_R3L3_abb R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(b!6) +#Molecules P_T2_bL_R3L3_baa R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6) +#Molecules P_T2_bL_R3L3_bab R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(b!6) +#Molecules P_T2_bL_R3L3_bba R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6) +#Molecules P_T2_bL_R3L3_bbb R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(b!6) + +Molecules T2aS1 R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4) +Molecules T2aS2 R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4) +Molecules T2aS3 R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4) +Molecules T2aS4 R2(a!+,b!1,aFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4) + +Molecules T2aL1 R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL2 R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL3 R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL4 R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL5 R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL6 R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL7 R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) +Molecules T2aL8 R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6,b!+,aFlag~1) + +Molecules T2bS1 R2(a!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +Molecules T2bS2 R2(b!1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +Molecules T2bS3 R2(a!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) +Molecules T2bS4 R2(b!1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!+,bFlag~1) + +Molecules T2bL1 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6) +Molecules T2bL2 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(b!6) +Molecules T2bL3 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6) +Molecules T2bL4 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(a!2,b!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(b!6) +Molecules T2bL5 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(a!6) +Molecules T2bL6 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(a!4,b!5).L2(r!5,r!6).R2(b!6) +Molecules T2bL7 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(a!6) +Molecules T2bL8 R2(a!+,b!1,bFlag~1).L2(r!1,r!2).R2(b!2,a!3).L2(r!3,r!4).R2(b!4,a!5).L2(r!5,r!6).R2(b!6) + +######## +# T3=3 # +######## + + + +end observables + +begin functions + +# global functions + +#avg_agg_size_R0()=(R0tot-monomericR0)/(R0_naggs+epsilon) +avg_agg_size_R1()=(R1tot-monomericR1)/(R1_naggs+epsilon) +avg_agg_size_R2()=(R2tot-monomericR2)/(R2_naggs+epsilon) +#avg_agg_size_R3()=(R1tot-monomericR3)/(R3_naggs+epsilon) + +# local functions + +######## +# T2=2 # +######## + +mobility_factor_T2(x) if(R2tot(x)>T2,0,1) +large_T2(x) if(R2tot(x)>T2,1,0) + +twobigprod_T2a(x) if((T2aS1(x)>0 || T2aS2(x)>0 || T2aS3(x)>0 || T2aS4(x)>0) && \ +(T2aL1(x)>0 || T2aL2(x)>0 || T2aL3(x)>0 || T2aL4(x) >0 || T2aL5(x)>0 || T2aL6(x)>0 || T2aL7(x)>0 || T2aL8(x)>0),1,0) + +twobigprod_T2b(x) if((T2bS1(x)>0 || T2bS2(x)>0 || T2bS3(x)>0 || T2bS4(x)>0) && \ +(T2bL1(x)>0 || T2bL2(x)>0 || T2bL3(x)>0 || T2bL4(x) >0 || T2bL5(x)>0 || T2bL6(x)>0 || T2bL7(x)>0 || T2bL8(x)>0),1,0) + +######## +# T3=3 # +######## + +#mobility_factor_T3(x) if(R3tot(x)>T3,0,1) +#large_T3(x) if(R3tot(x)>T3,1,0) + + + +end functions + +begin reaction rules + +########### +# Model 0 # +########### + +#L0(r,r)+R0(l)<->L0(r,r!1).R0(l!1) kf,kr +#L0(r,r!+)+R0(l)<->L0(r!1,r!+).R0(l!1) kpx,kmx + +########### +# Model 1 # - no immobilization +########### + +# capture/release of free ligand +L1(r,r)+R1(l1)<->L1(r,r!1).R1(l1!1) kf,kr +L1(r,r)+R1(l2)<->L1(r,r!1).R1(l2!1) kf,kr + +# ligand-mediated receptor crosslinking +L1(r,r!+)+R1(l1)<->L1(r!1,r!+).R1(l1!1) kpx,kmx +L1(r,r!+)+R1(l2)<->L1(r!1,r!+).R1(l2!1) kpx,kmx + +########### +# Model 2 # - immobilization of aggregates containing more than 2 receptors +########### + +# capture/release of free ligand +L2(r,r)+R2(a)<->L2(r!1,r).R2(a!1) kf,kr +L2(r,r)+R2(b)<->L2(r,r!1).R2(b!1) kf,kr + +# receptor crosslinking by tethered ligand with consideration of receptor aggregate mobility +# reactant x is mobile or immobile; reactant y is mobile +L2(r,r!+)+%y:R2(a)->L2(r!1,r!+).R2(a!1) kpx*mobility_factor_T2(y) +L2(r!+,r)+%y:R2(b)->L2(r!+,r!1).R2(b!1) kpx*mobility_factor_T2(y) +# reactant x is mobile; reactant y is immobile +%x:L2(r,r!+)+%y:R2(a)->L2(r!1,r!+).R2(a!1) \ +FunctionProduct("kpx*mobility_factor_T2(x)","1.0-mobility_factor_T2(y)") +%x:L2(r!+,r)+%y:R2(b)->L2(r!+,r!1).R2(b!1) \ +FunctionProduct("kpx*mobility_factor_T2(x)","1.0-mobility_factor_T2(y)") +# reactant x is immobible; reactant y is immobile +%x:L2(r,r!+)+%y:R2(a)->L2(r!1,r!+).R2(a!1) \ +FunctionProduct("sigma*kpx*large_T2(x)","large_T2(y)") +%x:L2(r!+,r)+%y:R2(b)->L2(r!+,r!1).R2(b!1) \ +FunctionProduct("sigma*kpx*large_T2(x)","large_T2(y)") + +L2(r!1,r!+).R2(a!1,aFlag~0)->L2(r!1,r!+).R2(a!1,aFlag~1) kmx +L2(r!+,r!1).R2(b!1,bFlag~0)->L2(r!+,r!1).R2(b!1,bFlag~1) kmx + +L2(r!1,r!+).R2(a!1,aFlag~1)->L2(r,r!+)+R2(a,aFlag~0) delta +L2(r!+,r!1).R2(b!1,bFlag~1)->L2(r!+,r)+R2(b,bFlag~0) delta + +# THE LINES BELOW DO NOT WORK BECAUSE OF A BUG IN BioNetGen/NFsim: +%z:L2(r!1,r!+).R2(a!1,aFlag~1)->L2(r!1,r!+).R2(a!1,aFlag~0) rho*twobigprod_T2a(z) +%z:L2(r!+,r!1).R2(b!1,bFlag~1)->L2(r!+,r!1).R2(b!1,bFlag~0) rho*twobigprod_T2b(z) + +########### +# Model 3 # +########### + + + +end reaction rules + +end model + +begin actions + +simulate({method=>"nf",gml=>2000000,complex=>1,get_final_state=>0,print_functions=>1,\ + t_start=>0,t_end=>300,n_steps=>300}) + +end actions diff --git a/tests/v1/test_bngl_corpus.py b/tests/v1/test_bngl_corpus.py new file mode 100644 index 00000000..2d479452 --- /dev/null +++ b/tests/v1/test_bngl_corpus.py @@ -0,0 +1,177 @@ +"""Corpus regression gate: the BNGL reader agrees with BNG2.pl on real models. + +Asserts :func:`petab.v1.models.bngl_model.parse_bngl` enumerates the same model +entities BNG2.pl does, over a curated set of public community BNGL models under +``bngl_corpus/``. BNG2.pl's answers are cached in ``bngl_corpus/golden.json``: +the entity name sets BNG2.pl emits from ``writeModel`` (its canonical parse, no +network generation). So the test itself needs **no BNG2.pl** -- it compares the +reader's output against the frozen oracle. Seed species are compared by +molecule composition, which absorbs BNG2.pl's pattern canonicalization (``t`` +vs ``t()``, component reordering, ``@compartment`` prefix vs suffix) while +still catching a genuinely missing or invented species. + +Regenerate the golden after adding a fixture or changing the reader (needs a +BNG2.pl on ``BNGPATH``/``PATH``) and review the diff:: + + python tests/v1/test_bngl_corpus.py +""" + +from __future__ import annotations + +import json +import re +import shutil +import subprocess +import sys +import tempfile +from pathlib import Path + +import pytest + +from petab.v1.models.bngl_model import parse_bngl + +_CORPUS = Path(__file__).parent / "bngl_corpus" +_GOLDEN = _CORPUS / "golden.json" +_MODELS = sorted(_CORPUS.glob("*.bngl")) + + +def _molecule_name(part): + """The molecule name in one ``.``-separated piece of a species pattern, + ignoring an ``@compartment:`` prefix, a ``$`` clamp, and any components. + + BNG2.pl writes the compartment prefix before the clamp (``@C::$ADP()``), so + the prefix is stripped first, then the ``$``. + """ + part = part.strip() + prefix = re.match(r"@\w+::?", part) + if prefix: + part = part[prefix.end() :] + part = part.lstrip("$") + name = re.match(r"(\w+)", part) + return name.group(1) if name else None + + +def seed_composition(patterns): + """Reordering-robust seed-species signature: each species as its sorted + molecule-name list, the whole sorted. JSON-safe (lists, not tuples).""" + sig = [] + for pattern in patterns: + mols = sorted( + n for n in (_molecule_name(p) for p in pattern.split(".")) if n + ) + sig.append(mols) + return sorted(sig) + + +def _reader_entities(entities): + """The entity name sets the golden records / the reader is checked against. + BNG2.pl-generated ``_``-names (a valid BNGL name starts with a letter) are + excluded from both sides.""" + + def named(names): + return sorted(n for n in names if not n.startswith("_")) + + return { + "parameters": named(entities.parameters), + "observables": sorted(entities.observable_names), + "functions": named(entities.function_names), + "molecule_types": sorted(entities.molecule_type_names), + "compartments": sorted(entities.compartment_names), + "seed_composition": seed_composition(entities.seed_species), + } + + +@pytest.mark.parametrize("model", _MODELS, ids=lambda p: p.stem) +def test_reader_matches_bng2_golden(model): + golden = json.loads(_GOLDEN.read_text()) + expected = golden[model.name] + actual = _reader_entities(parse_bngl(model.read_text())) + assert actual == expected, ( + f"{model.name}: reader disagrees with the BNG2.pl golden -- " + f"regenerate with `python {__file__}` if the change is intended" + ) + + +# -------------------------------------------------------------------------- +# Golden regeneration -- dev-only, needs BNG2.pl. Runs `writeModel` on each +# model's definition blocks (actions stripped, so no network generation), +# re-parses the canonical BNGL BNG2.pl emits, and records its entity name sets. +# -------------------------------------------------------------------------- + +_MODEL_BLOCKS = frozenset( + { + "parameters", + "molecule types", + "molecules", + "seed species", + "species", + "observables", + "functions", + "compartments", + "reaction rules", + "rules", + "energy patterns", + "population types", + "population maps", + } +) + + +def _strip_to_model(text): + """The model-definition blocks plus a single ``writeModel`` action; drops + actions (even nested in ``begin model``) and bare top-level directives.""" + out, stack = [], [] + for line in text.splitlines(): + s = line.split("#", 1)[0].strip() + begin = re.match(r"begin\s+(.+)", s, re.I) + end = re.match(r"end\s+(.+)", s, re.I) + if begin: + stack.append(begin.group(1).strip().lower()) + if stack[-1] != "actions": + out.append(line) + elif end: + top = stack.pop() if stack else None + if top != "actions": + out.append(line) + elif "actions" not in stack and [b for b in stack if b != "model"]: + out.append(line) + return "\n".join(out) + '\nwriteModel({prefix=>"canon"})\n' + + +def _canonical_bngl(model_text, bng2): + work = tempfile.mkdtemp(prefix="bnggold_") + try: + (Path(work) / "in.bngl").write_text(_strip_to_model(model_text)) + result = subprocess.run( # noqa: S603 + [bng2, "in.bngl"], + cwd=work, + capture_output=True, + text=True, + timeout=120, + ) + canon = Path(work) / "canon.bngl" + if result.returncode != 0 or not canon.is_file(): + raise RuntimeError( + f"BNG2.pl rejected the model:\n{result.stdout}{result.stderr}" + ) + return canon.read_text(encoding="utf-8", errors="replace") + finally: + shutil.rmtree(work, ignore_errors=True) + + +def _regenerate(): + from petab.v1.models.bngl_model import _locate_bng2 + + bng2 = _locate_bng2() + if bng2 is None: + sys.exit("No BNG2.pl found (set BNGPATH or put it on PATH).") + golden = {} + for model in _MODELS: + canon = _canonical_bngl(model.read_text(), bng2) + golden[model.name] = _reader_entities(parse_bngl(canon)) + _GOLDEN.write_text(json.dumps(golden, indent=1, sort_keys=True) + "\n") + print(f"wrote {_GOLDEN} ({len(golden)} models)") + + +if __name__ == "__main__": + _regenerate() diff --git a/tests/v1/test_model_bngl.py b/tests/v1/test_model_bngl.py new file mode 100644 index 00000000..b4641faf --- /dev/null +++ b/tests/v1/test_model_bngl.py @@ -0,0 +1,340 @@ +"""Tests for petab.v1.models.bngl_model""" + +from pathlib import Path + +import pytest + +from petab.v1.models.bngl_model import BnglEntities, BnglModel, parse_bngl + +#: A self-contained BNGL model exercising every entity kind read by the +#: PEtab layer: parameters, a molecule type, a seed species, an observable, +#: and a global function (the parabola measurement model). +EXAMPLE_BNGL = Path(__file__).parent / "bngl" / "parabola.bngl" + + +@pytest.fixture +def model() -> BnglModel: + return BnglModel.from_file(EXAMPLE_BNGL) + + +def test_parse_bngl(): + entities = parse_bngl(EXAMPLE_BNGL.read_text()) + assert isinstance(entities, BnglEntities) + assert entities.parameters == { + "v1": "5", + "v2": "5", + "v3": "5", + "k_cond": "1", + } + assert entities.observable_names == {"x"} + assert entities.function_names == {"y"} + assert entities.molecule_type_names == {"counter"} + assert entities.seed_species == {"counter()"} + assert entities.compartment_names == frozenset() + + +def test_parameter_ids_and_values(model): + assert set(model.get_parameter_ids()) == {"v1", "v2", "v3", "k_cond"} + assert model.get_parameter_value("v1") == 5.0 + assert dict(model.get_free_parameter_ids_with_values()) == { + "v1": 5.0, + "v2": 5.0, + "v3": 5.0, + "k_cond": 1.0, + } + + +def test_get_parameter_value_unknown_raises(model): + with pytest.raises(ValueError): + model.get_parameter_value("nope") + + +def test_expression_valued_parameter_is_not_evaluated(): + # A numeric RHS coerces to float; an expression RHS is confined to + # NotImplementedError rather than evaluated (that needs BNG2.pl). + entities = parse_bngl( + "begin parameters\n base 2\n k_on 2*base\nend parameters\n" + ) + model = BnglModel(entities, model_id="m") + assert model.get_parameter_value("base") == 2.0 + with pytest.raises(NotImplementedError): + model.get_parameter_value("k_on") + # The expression-valued parameter is still an enumerated entity, but it + # contributes no introspection-grade value. + assert dict(model.get_free_parameter_ids_with_values()) == {"base": 2.0} + + +# -- grammar hardening: block aliases + seed-species "$" clamp --------------- +# Kept in sync with PyBNF's sibling reader (pybnf/petab/_bngl.py, ADR-0026); +# these cases are the anchor that keeps the two block scanners from drifting. + + +def test_seed_species_dollar_clamp_is_stripped(): + # SeedSpeciesDefn = ["$"], Species, WS, MathExpression -- the "$" fixes the + # concentration but is not part of the species identity, so the enumerated + # state variable is the bare pattern (attached "$A()" and spaced "$ A()"). + entities = parse_bngl( + "begin seed species\n $A() 100\n $ B() 0\n C() 5\nend seed species\n" + ) + assert entities.seed_species == frozenset({"A()", "B()", "C()"}) + assert "$A()" not in entities.seed_species # the marker never leaks + + +def test_rejected_block_aliases_are_not_honored(): + # The grammar doc lists `molecules`/`rules` as aliases, but BNG2.pl 2.9.3 + # REJECTS both ("Could not process block type"); only `species` is real. To + # match the reference implementation the reader must NOT treat + # `begin molecules`/`begin rules` as their canonical blocks. + entities = parse_bngl("begin molecules\n A()\n B(x)\nend molecules\n") + assert entities.molecule_type_names == frozenset() + + +def test_seed_species_block_alias(): + # `begin species` is BNG's short alias for `begin seed species` -- and the + # "$" clamp is stripped under the alias spelling too. + entities = parse_bngl("begin species\n $A() 100\n B() 0\nend species\n") + assert entities.seed_species == frozenset({"A()", "B()"}) + + +def test_line_continuation_is_joined(): + # A trailing "\" continues the logical line (BNG2.pl readFile). Without + # joining, a continued parameter reads as the value "\"; the join + # concatenates directly -- no space -- so a token split across the break + # ("1e\"+"3" -> "1e3") rejoins, matching BNG2.pl. + entities = parse_bngl( + "begin parameters\n" + " minusb = \\\n(p4-1)/(p4*(1+p2))\n" # continued expression value + " r 1e\\\n3\n" # token split -> 1e3, no space + " a = 1+\\\n2+\\\n3\n" # chained continuation + "end parameters\n" + ) + assert entities.parameters["minusb"] == "(p4-1)/(p4*(1+p2))" + assert entities.parameters["r"] == "1e3" + assert entities.parameters["a"] == "1+2+3" + + +def test_backslash_in_comment_is_not_a_continuation(): + # BNG2.pl strips the comment before testing for a trailing "\", so a "\" + # inside a comment must not swallow the next line. + entities = parse_bngl( + "begin parameters\n k 1 # note \\\n j 2\nend parameters\n" + ) + assert entities.parameters == {"k": "1", "j": "2"} + + +def test_indexed_declarations(): + # Legacy .net-style leading index (LineLabel = {Digit}, WS): the index must + # not be read as the name. + entities = parse_bngl( + "begin parameters\n 1 L0 1\n 2 R0 2\nend parameters\n" + "begin seed species\n 1 A() 100\n 2 B() 50\nend seed species\n" + ) + assert entities.parameters == {"L0": "1", "R0": "2"} + assert entities.seed_species == frozenset({"A()", "B()"}) + + +def test_labeled_seed_species(): + # Named line label (LineLabel = Name, ":"): "CD14: CD14(...)" -- the label, + # which here even equals the molecule name, must not be read as the + # species. The label is stripped before the "$" clamp. + entities = parse_bngl( + "begin seed species\n" + " CD14: CD14(TLR4,MD2) v1\n" + " clamp: $MD2(x~0) v2\n" + "end seed species\n" + ) + assert entities.seed_species == frozenset({"CD14(TLR4,MD2)", "MD2(x~0)"}) + + +def test_line_label_does_not_over_strip(): + # A normal `name value` param and an `@compartment:` species must be left + # alone -- a compartment prefix carries "@", so it is not a bare label. + entities = parse_bngl( + "begin parameters\n NA = 6.02e23\n k1 1.0\nend parameters\n" + "begin seed species\n @PM:Rec() 100\nend seed species\n" + ) + assert entities.parameters == {"NA": "6.02e23", "k1": "1.0"} + assert entities.seed_species == frozenset({"@PM:Rec()"}) + + +def test_alias_does_not_shadow_the_canonical_block(): + # The `species` alias must not swallow the block whose name it is a + # substring of: `seed species` and `molecule types` stay distinct. + entities = parse_bngl( + "begin molecule types\n Counter()\nend molecule types\n" + "begin seed species\n $Counter() 1\nend seed species\n" + ) + assert entities.molecule_type_names == frozenset({"Counter"}) + assert entities.seed_species == frozenset({"Counter()"}) + + +def test_state_variable_ignores_the_clamp(): + # A clamped seed species is still a state variable under its bare id + # (is_state_variable drives CheckModel's species cross-checks). + model = BnglModel( + parse_bngl("begin seed species\n $A() 100\nend seed species\n"), + model_id="m", + ) + assert model.is_state_variable("A()") + assert not model.is_state_variable("$A()") + + +def test_has_entity_spans_full_declared_namespace(model): + # parameter, observable, global function, molecule type, seed species. + for entity in ("v1", "x", "y", "counter", "counter()"): + assert model.has_entity_with_id(entity) + # Prefixed PEtab IDs and unknowns are not model entities. + for non_entity in ("obs_x", "func_y", "nope"): + assert not model.has_entity_with_id(non_entity) + + +def test_symbol_allowed_is_the_paramlist_only(model): + # parameters u observables u global functions (the BNG ParamList). + for symbol in ("x", "y", "v1"): + assert model.symbol_allowed_in_observable_formula(symbol) + # A molecule type / seed species is an entity but not a formula symbol. + for non_symbol in ("counter", "counter()", "nope"): + assert not model.symbol_allowed_in_observable_formula(non_symbol) + + +def test_is_state_variable_is_seed_species_only(model): + assert model.is_state_variable("counter()") + assert not model.is_state_variable("v1") + assert not model.is_state_variable("x") + + +def test_valid_ids_for_condition_table_is_params_and_compartments(): + entities = parse_bngl( + "begin parameters\n k 1\nend parameters\n" + "begin compartments\n EC 3 1.0\nend compartments\n" + ) + model = BnglModel(entities, model_id="m") + assert set(model.get_valid_ids_for_condition_table()) == {"k", "EC"} + assert model.has_entity_with_id("EC") + # A compartment is an entity but not an observable-formula symbol. + assert not model.symbol_allowed_in_observable_formula("EC") + + +def test_invalid_model_id_raises(): + with pytest.raises(ValueError, match="not a valid identifier"): + BnglModel(parse_bngl(""), model_id="1nope") + + +def test_repr(model): + assert repr(model) == "" + + +def test_is_valid(model): + # A valid model validates: a real ``BNG2.pl --check`` where a BNG backend + # is locatable, or a graceful ``True`` fallback where it is not. + assert model.is_valid() is True + # A buffer-loaded model has no file to check and falls back to True. + assert BnglModel(parse_bngl(""), model_id="m").is_valid() is True + + +def test_is_valid_detects_broken_model_when_bng_available(tmp_path): + from petab.v1.models.bngl_model import _locate_bng2 + + if _locate_bng2() is None: + pytest.skip("BNG2.pl not available; is_valid falls back to True") + # An undefined parameter in the function body -> BNG2.pl --check fails. + broken = tmp_path / "broken.bngl" + broken.write_text( + "begin model\n" + " begin parameters\n k 1\n end parameters\n" + " begin functions\n f()=k*undefined_symbol\n end functions\n" + "end model\n" + ) + assert BnglModel.from_file(broken).is_valid() is False + + +def test_model_factory_routes_bngl(): + from petab.v1.models.model import model_factory + + model = model_factory(EXAMPLE_BNGL, "bngl") + assert isinstance(model, BnglModel) + assert model.type_id == "bngl" + + +def test_full_petab_validation_is_clean(tmp_path): + """A ``language: bngl`` problem loads via ``Problem.from_yaml`` and passes + every default validation task -- including the model-cross checks that + read the BnglModel (CheckModel, CheckObservablesDoNotShadowModelEntities, + CheckAllParametersPresentInParameterTable, CheckValidConditionTargets, + CheckInitialChangeSymbols). + """ + petab_v2 = pytest.importorskip("petab.v2") + from petab.v2 import Problem + from petab.v2.lint import ( + ValidationIssueSeverity, + default_validation_tasks, + ) + + # The BNGL model is the only genuinely BNGL-specific artifact; the PEtab + # tables are built with the v2 API so their format is guaranteed correct. + (tmp_path / "parabola.bngl").write_text(EXAMPLE_BNGL.read_text()) + + problem = Problem() + for param in ("v1", "v2", "v3"): + problem.add_parameter(param, estimate=True, lb=0, ub=10) + # observableFormula is the bare model name (x is an observable, y a global + # function); both are vouched for by the model, so neither is demanded as + # an output parameter in the parameter table. + problem.add_observable( + "obs_x", "x", noise_formula="0.1", noise_distribution="normal" + ) + problem.add_observable( + "func_y", "y", noise_formula="0.1", noise_distribution="normal" + ) + # A condition targets a model parameter that is *not* estimated -- a valid + # condition target (parameters u compartments) that exercises + # CheckValidConditionTargets against the BnglModel. + problem.add_condition("cond1", k_cond=2) + problem.add_experiment("exp1", 0, "cond1") + for t in (0.0, 1.0, 2.0): + problem.add_measurement( + "obs_x", time=t, measurement=1.0, experiment_id="exp1" + ) + problem.add_measurement( + "func_y", time=t, measurement=2.0, experiment_id="exp1" + ) + + petab_v2.write_parameter_df( + problem.parameter_df, tmp_path / "parameters.tsv" + ) + petab_v2.write_observable_df( + problem.observable_df, tmp_path / "observables.tsv" + ) + petab_v2.write_condition_df( + problem.condition_df, tmp_path / "conditions.tsv" + ) + petab_v2.write_experiment_df( + problem.experiment_df, tmp_path / "experiments.tsv" + ) + petab_v2.write_measurement_df( + problem.measurement_df, tmp_path / "measurements.tsv" + ) + + (tmp_path / "problem.yaml").write_text( + "format_version: 2.0.0\n" + "parameter_files: [parameters.tsv]\n" + "model_files:\n" + " parabola:\n" + " location: parabola.bngl\n" + " language: bngl\n" + "condition_files: [conditions.tsv]\n" + "experiment_files: [experiments.tsv]\n" + "observable_files: [observables.tsv]\n" + "measurement_files: [measurements.tsv]\n" + ) + + loaded = Problem.from_yaml(str(tmp_path / "problem.yaml")) + assert isinstance(loaded.model, BnglModel) + + errors = [ + (type(task).__name__, issue.message) + for task in default_validation_tasks + if (issue := task.run(loaded)) is not None + and issue.level == ValidationIssueSeverity.ERROR + ] + assert errors == []