Integrated model description

Model construction

BioUML platform for model reconstruction

BioUML (Biological Universal Modeling Language) (Kolpakov et al., 2019) [1] is an integrated platform for modeling of biological systems and developed in Java. The tool is applicable to solve a wide range of tasks including access to diverse biological databases, mathematical description and visual representation of biological systems, numerical calculations and parametric and other types of model analysis. The main features of the BioUML are:

• an ability to use both standalone version of the tool and web-services remotely;
• a support of generally accepted standards for model description of biological systems (SBML) and their graphical notation (SBGN) (Hucka et al., 2018; Le Novere et al., 2009)[2][3];
• visual modeling of biological systems and processes: an user has an opportunity to graphically construct and edit the developed model;
• a support of different mathematical representations (ordinary differential equations, algebraic equations, discrete events and stochastic modeling);
• a modular platform architecture that facilitates extension and/or addition of new types of models, methods of numerical calculations etc.

Visual modeling

Representation of investigating systems as graphical diagrams by means of a software supporting visual modeling can significantly facilitate the procedure of the model reconstruction. We consider visual modeling as a formal graphical representation of the system and/or modeling processes as a diagram and consequent dynamic simulation based on the representation. The graphical notation is a crucial component of visual modeling that allows one to formally and completely build a model. A visual model can be presented by some types of diagrams enabling description of diverse aspects of the structure and function of a complex system with different levels of details. This formal graphical representation is a basis for automatic code generation by specialized tools to simulate the model. We have made an extension of well-known SBGN notation (Le Novere, 2009)[3] in order to simulate physiological processes which require not only description between metabolites, but also ability to use algebraic, differential equations and instant transition of the system from one state to another. In addition, connections between equations indicate signal transduction in the model while interface ports of modules (or submodels) have also a direction (input, output or contact). A meta-model is a basis of the visual modeling in BioUML which ensures a formalism for complex description, graphical representation and numerical simulation of biological systems on different levels of their hierarchical organization. A meta-model consists of three interrelated levels of complex systems description:

• graphical representation - system’s structure is described as a compartmentalized graph;
• database - each element of the graph can include a reference on a certain object in the database;
• runned model - an element of the model (variables, mathematical equations, discrete events, states and transitions) can be associated with an element of the graph (vertices, arcs and compartments). As an example, vertices of the graph can be represented by variables or states of the system, while arcs of the graph correspond to equations describing changes of these variables or transitions between two states.

The description of the biological system as a meta-model is used to generate a Java code reflecting the model as a system of algebraic and/or differential equations, considering delay components, piecewise functions, discrete events and transitions. To generate a code the specific simulation engine is employed which defines the model type and corresponding simulation method.

A multi-compartmental complex model

A BioUML diagram describing a modular multi-compartmental model contains interconnected elements or modules (submodels) each of which is referred to another diagram (also may be modular). A directed connection between input and output nodes determines the signal transduction from a module to another one, while undirected relation between contacts reflects signals exchange between modules (Figure 1). According to this methodology, an integrated mathematical model describing energy metabolism of the human skeletal muscle (Kiselev et al., 2019, Akberdin et al., 2020)[4][5] considering Ca²⁺-dependent signaling pathway and downstream regulatory processes of early and late response genes expression has been built. A complex mathematical model (Akberdin et al., 2020)[5] developed by Li and coathours in MATLAB has been rebuilt in BioUML as an initial model of the energy metabolism of the human skeletal muscle taking into account quantitative differences between fiber I and II types (Figure 2). An activation mechanism that enhances energy metabolism via transport and reaction fluxes due to physical exercise was harnessed as the stress function depending on general work rate parameter. The work rate parameter defines intensity of the physical exercise and variates depending on the mode of the exercise.

Figure 1. An integrated modular model linking metabolism, Ca²⁺-dependent signaling transduction and regulation of gene expression in human skeletal muscle. Grey arrows in the metabolic module represent transport reactions of metabolites, blue arrows mean external control function like consideration of physical exercise. Module «Muscle tissue» consists of submodules «Cytosol» and «Mitochondrion», which in turn contain equations describing enzymatic reactions inside the certain compartment (А). A certain mathematical module has been developed for each early and late response gene (B).

Upgrade of the model (metabolic level)

It is worth to note that values of activation coefficients associated with ATPase (Stienen et al., 1996, He et al., 2000, Szentesi et al., 2001, Barclay, 2017)[6][7][8][9] and pyruvate dehydrogenase reaction fluxes for type I and type II fibers (Parolin et al., 1999, Kiilerich et al., 2008, Albers et al., 2015)[10][11][12] as well as time constant of ATPase flux rate coefficient in response to exercise were modified according to recently published data and estimations (Broxterman et al., 2017, Bartlett et al., 2020)[13][14]. Despite overall net glycogen breakdowns during muscle contraction, exercise also increases the activity of glycogen synthase (GS) (Wojtaszewski et al., 2001, Nielsen et al., 2003, Jensen et al., 2009, Jensen et al., 2012)[15][16][17][18]. The GS reaction results in ATP consumption, therefore GS reaction fluxes were modified according to (Wojtaszewski et al., 2001, Jensen et al., 2009, Jensen et al., 2012b)[15][17][18]. The rates of muscle glycogen synthesis during exercise assumed to be equal in type I and type II fibres and were estimated from average post-exercise glycogen synthesis data (Casey et al., 1995)[19]. To consider the allosteric regulation of AMPK activity (in corresponding modules) concentrations of free ADP and AMP in the cytosol were calculated using intracellular Cr, PCr, ATP and H+ concentrations as well as the equilibrium constants for creatine phosphokinase and adenylate kinases in each fiber type as described previously (Lawson et al., 1979, Dudley et al., 1985, Mannion et al., 1993)[20][21][22] (Figure 3).

Figure 2. A general SBGN diagram of the modular model describing metabolism in human muscle fibers (A) taking into account metabolic processes in the cytoplasm (B), in the mitochondrion (C) and transport reactions between two compartments (D).

The diagram of the modular model describing the metabolism of human skeletal muscle is presented in Figure 2. The cytosol includes metabolic reactions of the glycolysis, glycogenolysis and lipids metabolism, while tricarboxylic acid (TCA) cycle, ß-oxidation and oxidative phosphorylation reactions are presented in the mitochondria. The intermediate compartment between those is a transport module that contains passive and facilitated transport reactions for model intracellular species. Kinetic laws presenting metabolic and transport flux expressions exactly match the initial model developed by Li and coathors [23].

Signaling level

The concentration of Ca²⁺ ions in the myoplasm increases in proportion to the intensity of exercise. Ca²⁺ binds to calmodulin, thereby activating CaMKs and phosphatase calcineurin (Gehlert et al., 2015)[24]. CaMKII is the most abundant isoform in the human skeletal muscle, whereas CaMKI and CaMKIV are not expressed at detectable levels (Rose et al., 2006)[25]. An increase in CaMKII activity results in CREB1 Ser133 phosphorylation leading to activation of the transcription factor (Johannessen&Moens, 2007, Olesen et al., 2010) [26][27]. Calcineurin can dephosphorylate (and activate) CRTCs at Ser171 (CREB-regulated transcription coactivators) playing a key role in regulating the transcriptional activity of CREB1 (Altarejos et al., 2011)[28]. Another target of calmodulin is calcium/calmodulin-dependent protein kinase kinase 2 (CAMKK2) that phosphorylates AMPK Thr172 thereby activating the kinase (Abbott et al., 2009)[29]. In turn, activated AMPK can phosphorylate CREB1 Ser133 (Thomson et al., 2008)[30]. Collectively, these findings drove us to include in our model the Ca²⁺-dependent regulation of calmodulin, CREB1 (via CaMKII), CRTC (via calcineurin), and AMPK (via CaMKK2) (Figure 3). The amount of these proteins in human skeletal muscle was estimated using published proteomics and transcriptomics data (Murgia et al., 2017, Popov et al., 2019)[31][32] (see Supplementary data in Akberdin et al., 2020[5]). There are three different heterotrimeric complexes in the human skeletal muscles: α2β2γ1, α2β2γ3, and α1β2γ1 (Wojtaszewski et al., 2005)[33]. Distinct kinetic properties (an intrinsic enzyme activity, binding affinities of AMP, ADP and ATP to the specific isoform, sensitivity to de- and phosphorylation of AMPK heterotrimers) (Rajamohan et al., 2016, Ross et al., 2016)[34][35] and their subcellular localization (Pinter et al., 2013)[36] cause a differential regulation of the AMPK heterotrimers in vivo. The α2β2γ3 complex is phosphorylated and activated during moderate- to high-intensity exercise, while the activity associated with the other two AMPK heterotrimers is almost unchanged (Birk et al., 2006)[37]. However, the basal activity of α2β2γ3 complex is significantly lower than others. Taking into account the general AMPK basal and exercise-induced activity is considered as a sum of isoforms activities, all isoforms in the corresponding module was considered to quantitatively fit an experimental data obtained at baseline and after an exercise (Birk et al., 2006, Willows et al., 2017)[37][38]. AMPK is regulated by various ways: an up-stream kinase LKB1 can phosphorylate AMPK at Thr172 (Lizcano et al., 2004, Jansen et al., 2009)[39][40]. On the other hand, both ATP and AMP allosterically regulate AMPK: an exercise-induced decrease in intramuscular ATP increases its activity, while an increase in AMP activates it (Hardie et al., 2016, Li et al., 2017)[41][42]. Hence, in our model the AMPK is regulated via AMP, ATP, and LKB1, as well as CaMKK2 (as mentioned above).

Figure 3. A SBGN diagram of the Ca²⁺- and AMPK-dependent signaling pathway activating by contractile activity (aerobic exercise).

Gene expression level

An aerobic exercise induces expression of several hundreds of genes regulating many cell functions: energy metabolism, transport of various substances, angiogenesis, mitochondrial biogenesis. Regulation of the transcriptomic response to acute exercise includes dozens of transcription regulators [Popov et al. 2019] and seems to be extremely complex. Therefore to consider the response on gene expression level we exemplify the regulation of some genes encoding a transcription co-activator PGC-1α (encoded by PPARGC1A gene) and nuclear receptors NR4As - key exercise-induced regulators of the angiogenesis, mitochondrial biogenesis, fat and carbohydrate metabolism in skeletal muscle [Lira 2010, Pearen&Muscat 2018]. Expression of NR4A2, NR4A3 mRNA rapidly increases during the first hour after an aerobic exercise (early response genes) due activation of Ca2+\calcineurin-dependent signaling [8, 19, 46–48, Pearen&Muscat 2018]. We included in our model the Ca2+-dependent regulation (Ca2+\calcineurin-CaMKII-CREB1) of NR4As genes using data of contractile activity-specific mRNA response of these genes [8]. Expression of PPARGC1A mRNA rises 3 to 4 h after an exercise (late response gene МОЖЕТ БЫТЬ ЛУЧШЕ ВЕЗДЕ (и на рисунках) ЗАМЕНИТЬ НА gene with delay response Т.К. 4 ЧАСА СЛОЖНО НАЗВАТЬ ПОЗДНИМ ОТВЕТОМ) [8]. The transcription regulation of PPARGC1A via the canonical (proximal) and inducible (distal) promoters is very complicated, and includes Ca2+- and AMPK-dependent signaling, as well as CREB1 and its co-activator CRTC [Popov et al. 2015, Popov 2018]. The phosphorylation level of many signaling kinases drops to basal levels within the first hour after an aerobic exercise. Moreover, in a genome-wide study on various human tissues, it was shown that the phosphorylation level of CREB Ser133 does not always correlate with its transcriptional activity [53]. Therefore, we suggested the expression of late response genes (including PPARGC1A) is regulated by increasing the expression of one of the early response genes encoding transcription factors leading to a rapid increase corresponding protein (see Fig. 5 in Akberdin 2020). Analysis of contractile activity-specific transcriptomic data [8] showed that a rapid increase in the expression of genes encoding various TFs is observed already in the first hour after an exercise. It turned out that the binding motifs of some TFs (CREB-like proteins, as well as proteins of the AP-1 family: FOS and JUN) are located and intersected with each other both in the alternative and in the canonical promoters of the PPARGC1A gene [Akberdin 2020], i.e. these TFs can act as a potential regulators of this gene. This is consistent with the fact that these TFs can bind to DNA and regulate the expression of target genes as homo- and heterodimers [49, 50]. Based on these considerations, we included in the model the regulation of gene expression of early (NR4A2, NR4A3) and delayed (PPARGC1A) genes: early response genes are regulated via the activation of existent TFs (e.g. CREB1) and their co-activators (e.g. CRTC), while delayed response genes - via an increase in the expression of early response genes encoding transcription factors (transcription factor X in our model, Fig. N).

References

1. Kolpakov F, Akberdin I, Kashapov T, Kiselev L, Kolmykov S, Kondrakhin Y, Kutumova E, Mandrik N, Pintus S, Ryabova A, Sharipov R, Yevshin I, and Kel A. BioUML: an integrated environment for systems biology and collaborative analysis of biomedical data. Nucleic Acids Res. 2019 Jul 2;47(W1):W225-W233. DOI:10.1093/nar/gkz440 |
2. Hucka M, Bergmann FT, Dräger A, Hoops S, Keating SM, Le Novère N, Myers CJ, Olivier BG, Sahle S, Schaff JC, Smith LP, Waltemath D, and Wilkinson DJ. The Systems Biology Markup Language (SBML): Language Specification for Level 3 Version 2 Core. J Integr Bioinform. 2018 Mar 9;15(1). DOI:10.1515/jib-2017-0081 |
3. Le Novère N, Hucka M, Mi H, Moodie S, Schreiber F, Sorokin A, Demir E, Wegner K, Aladjem MI, Wimalaratne SM, Bergman FT, Gauges R, Ghazal P, Kawaji H, Li L, Matsuoka Y, Villéger A, Boyd SE, Calzone L, Courtot M, Dogrusoz U, Freeman TC, Funahashi A, Ghosh S, Jouraku A, Kim S, Kolpakov F, Luna A, Sahle S, Schmidt E, Watterson S, Wu G, Goryanin I, Kell DB, Sander C, Sauro H, Snoep JL, Kohn K, and Kitano H. The Systems Biology Graphical Notation. Nat Biotechnol. 2009 Aug;27(8):735-41. DOI:10.1038/nbt.1558 |
4. https://doi.org/10.17537/2019.14.373
5. https://doi.org/10.17537/2020.15.20