BioModelKit – An Integrative Framework for Multi-Scale Biomodel-Engineering

3.1 Module Terminology

In this section, we introduce the module terminology of the molecule-centred modularisation approach as in [32]. This formalisation is essential for the realisation of BMK including all of its features and functionalities. A brief example follows every term about the running example depicted in Figure 3. Figure 3(A) provides a short description of the molecular mechanisms, while Figure 3(B) depicts the corresponding modules, respectively the composed model.

A component is any chemical compound that contributes to a biological process at the molecular level.

A genetic component is a component with genetic information like genes, mRNAs, and proteins.

A non-genetic component is a component without genetic information like metabolites, second messengers, energy equivalents, ions, etc.

Figure 3:

(A) Molecular mechanism: Active kinase (K) phosphorylates its transcription factor (TF), the phosphorylated TF binds to the K gene, which starts the transcription, followed by protein translation of K. (B) Modular representation by PNs: Two protein modules for the TF and K, one gene, mRNA and protein degradation module for K according to the molecular mechanism in (A). Nodes highlighted in the same colour define a shared interface network. Modules share their tokens after model composition and thus, the dotted tokens appear on shared places enabling the execution of interactions.

A module M_c0 describes the functionality of one particular genetic component c₀, called main component, including its interactions with n_IC ≥ 0 components. Here, C(Mc0)={c0,…,cnIC} defines the total set of components of module M_c0. We distinguish mechanistic module types (protein, protein degradation, mRNA and gene modules) and causal module types (allelic influence and causal influence modules), see Section 3.2 for details.

A functional unit u is a defined part of a component c∈C(Mc0) with an assigned function that is of importance for the functionality of the component. A genetic component c∈C(Mc0) consists of n_FU ≥ 1 functional units. We assume that non-genetic components cannot be decomposed into more than one functional unit. Furthermore, U(c)={u1,…,unFU} defines the total set of functional units of component c∈C(Mc0), U(Mc0)=⋃c∈C(Mc0)U(c) defines the total set of functional units in module M_c0.

A molecular state s describes a specific physical constitution of a functional unit u∈U(Mc0). A functional unit u ∈ U(c) can adopt n_MS ≥ 1 different molecular states. Furthermore, S(u)={s1,…,snMS} defines the total set of molecular states of a functional unit u ∈ U(c) and S(Mc0)=⋃u∈U(Mc0)S(u) defines the summation of all sets S(u) over u∈U(Mc0). Each molecular state s∈S(Mc0) can be mapped to at least one functional unit according to the function λu(s)={u∈U(Mc0)∣s∈S(u)} and to at least one component according to the function λc(s)={c∈C(Mc0)∣s∈⋃u∈U(c)S(u)}.

An interaction state s_IS is a molecular state sIS∈S(Mc0) if it can be mapped to more than one component |λc(s)|>1 defining a complex k=λc(s). Accordingly, SIS(Mc0)={s∈S(Mc0)∣|λc(s)|>1}, SIS(Mc0)⊆S(Mc0) defines the total set of interaction states.

The term active molecular state defines the current state of a functional unit u∈U(Mc0), only one molecular state can be active at each time point. The function α:S(u)→{0,1} determines the active molecular state with the constraint ∑s∈S(u)α(s)=1. If α(s) = 1 the molecular state s is active, otherwise if α(s) = 0 the molecular state s is inactive.

The initial state s_0,u ∈ S(u) with α(s_0,u) = 1 and s_0,u ∉ SIS(Mc0) is assumed to be the ground state of a functional unit, which can not be an interaction state. This assumption is necessary to exclude inconsistencies during the composition of models. The union of all initial states is given by S0(Mc0)=⋃u∈U(c0)s0,u.

A molecular event e represents an action changing the molecular states of the functional units. The relations between molecular states in S(Mc0) are described by n_ME ≥ 1 molecular events. Here, E(Mc0)={e1,…,enME} defines the total set of molecular events.

The stoichiometric coefficientν(e,s)∈N0 (products) and ν(s,e)∈N0 (educts) defines the stoichiometry with which a molecular state s∈S(Mc0) mediates a molecular event e∈E(Mc0). The set of molecular states participating in a molecular event e can be distinguished into a set of educts ∙e={s∈S(Mc0)∣ν(s,e)≠0}, and a set products e∙={s∈S(Mc0)∣ν(e,s)≠0}. Vice versa, ∙s={e∈E(Mc0)∣ν(e,s)≠0} defines the set of molecular events, where the molecular state s is a product, and s∙={e∈E(Mc0)∣ν(s,e)≠0} defines the set of molecular events, where molecular state s is an educt. A molecular event e∈E(Mc0) can be defined as e={(ε(∙e),ε(e∙)),ε(∙e)→ε(e∙)}, where ε(∙e)=∑s∈∙eν(s,e)⋅s and ε(e∙)=∑s∈e∙ν(e,s)⋅s. Each molecular event e∈E(Mc0) occurs with a rate r(e).

An interaction event e_IS is a molecular event eIS∈E(Mc0) that involves more than one component, |λc(∙e∪e∙)|>1. The total set of interaction events is given by EIA(Mc0)={e:|λc(∙e∪e∙)|>1}.

A molecular process E(u) defines the set of molecular events changing the molecular state of a functional unit u ∈ U(c), such that E(u)=⋃u∈S(u)((∙s∪s∙)∖(∙s∩s∙)).

3.2 Module Types

The modularisation concept makes use of different module types to integrate molecular mechanisms and correlations on all of the main omic levels (metabolome, proteome, transcriptome, genome). Accordingly, we distinguish mechanistic module types [33], compare Figure 2 – Box 2:

A gene moduleMg,c0 represents the transcriptional activity of a gene controlled by the formation of the regulatory landscape and the pre-initiation complex. These processes include complex non-covalent interactions with transcription factors or other regulatory proteins interacting with a silencer or enhancer sequences of the gene.

A mRNA moduleMm,c0 represents the biosynthesis of a particular mRNA by transcription of the respective gene; the post-transcriptional modification of the mRNA including capping, (alternative) splicing and polyadenylation; the translation into the proteins encoded by the processed mRNA; and the degradation of the mRNA including its potential control through proteins or small interfering anti-sense RNA molecules. These processes include covalent modification and non-covalent interactions of the mRNA.

A protein moduleMp,c0 describes the functionality of a particular protein (single polypeptide chain), including changes in the protein conformation, non-covalent interactions with other components and covalent modifications that regulate the functionality. Thus, a protein module represents the formation and cleavage of covalent and non-covalent bonds, as well as conformational changes of the protein structure.

A protein degradation moduleMd,c0 represents the degradation of a protein by proteolysis in lysosomes, ubiquitin-dependent degradation by the proteasome, or degradation by digestive enzymes or by any other possible mechanism that may lead to the degradation or inactivation of the protein. The post-translational proteolytic processing is described in the corresponding protein module.

The module types introduced above exclusively rely on known molecular mechanisms, even though some mechanistic details may not necessarily be considered. To integrate causal relationships if molecular mechanisms unknown, we introduce two more module types [33]:

An allelic influence moduleMai,c0 represents the effects of alleles (mutated versions of a gene) on molecular processes. In contrast to gene modules, the described effects are causal influences, which might be directly or indirectly mediated by unknown processes.

A causal influence moduleMci,c0 describes the influence of arbitrary entities, others than alleles, on molecular processes.

3.3 Module Petri Net

The formal description of a module M, see Section 3.1 has to be translated into a PN N={P,T,F,f,v,m0}. Places represent molecular states and transitions represent molecular events. The function ϱp→s:P→S(M) (ϱt→e:T→E(M)) maps a place p ∈ P (transition t ∈ T) to a molecular state s∈S(M) (molecular event e∈E(M)), both mappings are bijective.

The PN N(Mc0)={P,T,F,f,v,m0} of a module M_c0 is defined as:

• Set of places P=⋃c∈C(Mc0)Pc, where

  1. Pc=⋃u∈U(c)Pu is the set of places representing a component c∈C(Mc0) and Pu={p:ϱp→s(p)∈S(u)} is the set of places representing a functional unit u∈U(Mc0);

• Set of transitions T=⋃c∈C(Mc0)Tc, where

  1. Tc=⋃u∈U(c)Tu is the set of transitions related to a component c∈C(Mc0) and Tu={t:ϱt→e(t)∈E(u)} is the set of transitions related to a functional unit u∈U(Mc0);

• Set of arcs F:=FSA⊆(P×T)∪(T×P);

• Arc-weights f: F→N0, where

  1. ∀s∈S(Mc0),e∈E(Mc0) with s∈∙e: fSA(ϱp→s−1(s),ϱt→e−1(e))=ν(s,e) (input arc), and

  2. ∀s∈S(Mc0),e∈E(Mc0) with s∈e∙: fSA(ϱt→e−1(e),ϱp→s−1(s))=ν(e,s) (output arc);

• Set of firing rates v:T→H with H=⋃t∈Th(t), where:

  1. ∀e∈E(Mc0):h(ϱt→e−1(e))=r(e);

• Initial marking m₀: P→N0:

  1. M_c0 is a gene, protein, causal or allelic influence module:

     1. ∀s∈S0(Mc0):m0(ϱp→s−1(s))=nc0 (n_c0 is the assumed initial number of copies for component c₀, here we assume nc0=1), and

     2. ∀s∈S(Mc0)∖S0(Mc0):m0(ϱp→s−1(s))=0,

  2. M_c0 is a protein degradation or mRNA module:

     1. ∀s∈S(Mc0):m0(ϱp→s−1(s))=0.

The protein degradation module Md,c0 is a special case with respect to the use of arcs. Assuming a transition t ∈ T represents the molecular event of protein degradation e=ϱt→e(t) of the protein defined by module Mp,c0′, then the set of places representing:

• none interaction states: PMp,c0′nonIS=ϱp→s−1(⋃u∈U(c0′)S(u)∖SIS(Mp,c0′)) are connected with transition t using:

  1. reset arcs fXA(PMp,c0′nonIS,t)=1, and

  2. marking-dependent standard arcs fSA(t,PMp,c0′nonIS)=m0(PMp,c0′nonIS);

• interaction states: PMp,c0′IS=ϱp→s−1(⋃u∈U(c0′)S(u)∩SIS(Mp,c0′)) are connected with transition t using:

  1. inhibitory arcs fIA(PMp,c0′IS,t)=1.

This transformation ensures that proteins are only degraded if they are not interacting with other components and that only one copy is degraded.

3.4 Interface Networks

Interface networks are shared subgraphs of modules describing an interaction between the represented components or a functional relationship (e.g. links between translation, transcription or protein degradation), in the most trivial case interface networks consist only of a single place. Identical interface networks are used to couple modules and are therefore crucial for the model composition.

A set of modules I={M1,…,Mn} of components involved in one particular interaction share the interface network NI(M)={PI,TI,FI,fI,vI,m0I(M)}, where NI(M)⊆N(M), M ∈ I. All nodes in an interface network NI(M) are declared as logical (fusion) nodes. The set of places P^I and transitions T^I in an interface network 𝒩_I with e.g. I={Mc0,1,Mc0,2} depends on the type of interaction among the two modules Mc0,1 and Mc0,2, which can be categorised into:

• regulation – assuming Mc0,1 and Mc0,2 can be any module type:

  1. TI={t:ϱt→e(t)∈EIA(Mc0,1)∩EIA(Mc0,2)}, and

  2. PI={p:ϱp→s(p)∈∙(EIA(Mc0,1)∩EIA(Mc0,2))∪(EIA(Mc0,1)∩EIA(Mc0,2))∙};

• transcription – assuming Mc0,1 is a gene module and Mc0,2 is an mRNA module:

  1. TI=∅, and

  2. PI={p:ϱp→s(p)∈S(Mc0,1)∧ϱp→s(p) is a transcriptionally active state of c0,2};

• translation – assuming Mc0,1 is an mRNA module and Mc0,2 is a protein module:

  1. T^I = ∅, and

  2. PI={p:ϱp→s(p)∈S0(Mc0,2)};

• degradation – assuming Mc0,1 is a protein module and Mc0,2 is a protein degradation module:

  1. T^I = ∅, and

  2. PI⊆{p:ϱp→s(p)∈S(Mc0,1)}.

The redundant interface networks shared among modules might appear unnecessarily complicated, but they are of tremendous benefits for modules with complex interaction sites by securing the correct functioning of modules in a composed model. Even more, interface networks ensure that interactions can only be executed if all modules of components involved in the interaction are part of the composed model according to the real-world scenario, see next section.

5.1 Module Annotation

As discussed by Le Novère et al. [31] most published biological models are lost due missing access or insufficient characterisation. Based on these observations Le Novère et al. proposed a guideline for curating and encoding models called MIRIAM (Minimum Information Requested In the Annotation of Models). To adhere to the proposed MIRIAM guideline, a module consists not only of its underlying PN but also of an annotation file. We, therefore, defined the BMK markup language (BMKml), the XML Schema Definition (XSD) is documented in [32]. In summary, the annotation file specifies:

• module name, type and a cross-reference to the Ensembl database [34] of the corresponding component;

• information to clarify the authorship including contact details, creation and modification dates of the module, as well as a short description of the considered molecular mechanisms; and

• for each place, transition and parameter a description, literature references and cross-references to other databases giving more detail about the modelled molecular mechanisms or linking a GO annotation [35].

The module annotation provides a complete characterisation and documentation of the module to ease the understanding of the modelled molecular mechanisms and the reuse in composed models.

5.2 Module Construction

So far, we exploit three ways to obtain modules, compare Figure 2 – Box 1:

1. The direct engineering of modules in Snoopy [18] based on knowledge provided in textbooks, publications, and databases is the default approach to construct modules. We used this approach to compose models for the following case studies:

   • Pain signalling – the composed model comprises 38 modules of key players involved in molecular mechanisms of nociception, which occurs in the peripheral terminals of dorsal root ganglion neurones and is responsible for the detection of noxious and painful stimuli [36], [37], [38], [39], [40].

   • IL6-induced JAK-STAT signalling – the composed model comprises seven protein modules (JAK, STAT, IL6, IL6R, gp130, SHP2, SOCS3), as well as a protein degradation, mRNA and gene module for the feedback inhibitor SOCS3 [25], [41].

   • Phosphate regulatory network in enterobacteria – the composed model comprises nine protein modules (PhoA, PstS, PstA, PstB, PStC, PhoU, PhoR, PhoA), as well as a protein degradation, mRNA and gene module for PhoA, which sense and regulate the external Pi [33].

2. The reverse engineering of modules from complex (omic) datasets to obtain causal modules or map data to module prototypes [32], [33], e.g. as in case of conceptualised gene or mRNA modules, are straightforward approaches to integrate experimental results directly into composed models and thus, to test for their effect on the model’s behaviour.

3. The modularisation and thus, the integration of existing models will be another valuable resource to obtain modules. In particular, we suggest a routine to decompose mechanistic SBML encoded models [23] into a set of related modules. Modules can also be obtained from gene regulatory models defined by Boolean networks as described in [22].

5.3 Algorithmic Mutation

Life science showed in order to understand a biological system it is necessary to study not only the native system state but also its genetic variations. Genetic variations are a consequence of mutations in the genome, where some of the mutations are indifferent, and others lead to a gain or even a loss of function. Mutations result either naturally or are experimentally enforced. Modelling approaches need to introduce equivalent mutation workflows. BMK suggests algorithmic mutation based on the modularisation concept to conveniently mimic gene deletions and structure-function mutations, which are described in detail in [32].

Removing modules from a composed model mimics gene deletions. This mutation can be performed as single, double or even multiple gene knockouts to sufficiently block the effect of redundant gene functions.

Since functional units in a module have a relation to the structure of the modelled component, their alterations can mimic structure-function mutation. By deleting the transitions in the subgraph of a functional unit, the functional unit can only adopt a subset of the previously reachable molecular states and can therefore not execute its complete functionality. As stated above, such alteration could result in an indifferent effect, but it could also possibly increase or decrease a specific function encoded by the molecular mechanism.

Both mutation algorithms can be performed either on the entire set of modules or on a subset of modules/functional units that can be filtered through the provided references in the annotation file, see Section 5.1. The systematic mutation of a composed model yields a large set of alternative models and demands analysis workflows based on machine learning to reveal fundamental insights and identify structures with a significant effect on the system’s behaviour.

5.4 Spatiotemporal Transformation

Molecular mechanisms are affected by the spatial distribution and movement of the involved components. Therefore, the cellular arrangement concerning compartments, membranes and pools, as well as the cell geometry and size, have to be considered to understand the spatiotemporal behaviour of a system. Representing biomolecular mechanisms and movement in one coherent model is a challenging task.

Based on the modularisation concept, we propose a spatiotemporal transformation algorithm that integrates those spatial aspects into the composed model, see [24], [42] for details. The spatiotemporal transformation algorithm extends the composed model with a spatial model controlling the movement and interaction of the involved components without interfering the modelled molecular mechanisms. In the first step of the spatiotemporal transformation, the cellular space needs to be defined by a grid, as well as cellular substructures on that particular grid. This step includes the assignment of the components in the composed model to a set of substructures. The position of each component on the grid is encoded by a set of coordinate places (one, two or three places depending on the dimension of the grid) and the corresponding marking. The movement of the component is realised by increasing or decreasing the marking on the coordinate places by transitions with respect to the grid and its substructures. The introduced subnets responsible for the movement must further distinguish, whether the component is in an interactive state forming a complex or not, to guarantee a synchronous movement of all components forming an interactive complex. Additional rules have to be introduced to ensure that interaction between components encoded in the interface networks can only occur if they reside close to each other. The spatiotemporal transformation employs coloured PNs to allow for a scalable representation of the spatial model concerning the copy number of the involved components. The copy number of components has to be realised via colour and not via the marking, to distinguish the individual copies in respect to their position, movement and interaction.

Applying the spatiotemporal transformation to a composed model allows studying how space might affect the dynamic behaviour of a system without constituting a new additional model or framework only for this specific purpose.

Figure 4:

Basics of the spatial transformation concept.