A Modular Mathematical Model of Exercise-Induced Changes in Metabolism, Signaling, and Gene Expression in Human Skeletal Muscle

Skeletal muscle is the principal contributor to exercise-induced changes in human metabolism. Strikingly, although it has been demonstrated that a lot of metabolites accumulating in blood and human skeletal muscle during an exercise activate different signaling pathways and induce the expression of many genes in working muscle fibres, the systematic understanding of signaling–metabolic pathway interrelations with downstream genetic regulation in the skeletal muscle is still elusive. Herein, a physiologically based computational model of skeletal muscle comprising energy metabolism, Ca2+, and AMPK (AMP-dependent protein kinase) signaling pathways and the expression regulation of genes with early and delayed responses was developed based on a modular modeling approach and included 171 differential equations and more than 640 parameters. The integrated modular model validated on diverse including original experimental data and different exercise modes provides a comprehensive in silico platform in order to decipher and track cause–effect relationships between metabolic, signaling, and gene expression levels in skeletal muscle.


Introduction
Skeletal muscle tissue comprises about 40% of total body mass in lean adult humans and plays a crucial role in the control of whole-body metabolism and exercise tolerance. Regular low-intensity exercise (aerobic or endurance training) strongly increases vascular and mitochondrial density and oxidative capacity, improving fat and carbohydrate metabolism. These adaptations lead to an enhancement of muscle endurance performance and reduce the risk associated with the morbidity and premature mortality of chronic cardiovascular and metabolic diseases [1,2].
Acute aerobic exercise induces significant metabolic changes in the working skeletal muscle, which in turn activate numerous signaling molecules. Changes in the content of Ca 2+ ions in skeletal muscle play a fundamental role in the regulation of the activity of contractile proteins and enzymes involved in energy metabolism. In addition, a contractioninduced increase in the content of Ca 2+ ions in the myoplasm significantly affects the • SBML-Systems Biology Markup Language [22] serves for a formal description of mathematical models. BioUML supports all versions of SBML from l1v2 to the latest l3v2, including "comp" [23]. • SBGN-Systems Biology Graphical Notation [24] is used for visual description of model elements (complexes, compartments, molecule types, reactions, etc.). BioUML completely supports SBGN Process Description diagrams and uses them to visually represent SBML models. We also support the XML markup language SBGN-ML (https://github.com/sbgn/sbgn/wiki/SBGN_ML accessed on 8 June 2021), which facilitates the exchange of SBGN diagrams between tools. • Antimony-human-readable text format, which supports most of the SBML features [25]. BioUML automatically processes it into an SBML diagram in SBGN notation. BioUML supports import and export into antimony format.
However, these standards are not sufficient for some tasks. Thus, we suggest extension of the SBGN Process Description diagram type and Antimony format and demonstrate how they can improve the construction of complex biological models using visual modelling. These extensions supported by the BioUML platform will be described below.

Visual Modelling
Representation of investigated systems as graphical diagrams by means of software supporting visual modeling can significantly facilitate the procedures of the model reconstruction.
Following a paradigm of visual modelling, a user creates mathematical models as visual diagrams. Each component of the diagram corresponds to a particular mathematical object of the model (variable, reaction, metabolite, equation, etc.). Users may additionally edit those elements by changing their properties (i.e., initial value of a variable, kinetic law for the reaction, etc.). Based on this visual representation as well as on defined properties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Element Name Glyph Description
Equation law for the reaction, etc.). Based on this visual representation as well as on defined properties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Event
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".

Function
Function receives argument values, calculates, and returns result of the given expression.

Constraint
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.

Tabular data
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function.
Event law for the reaction, etc.). Based on this visual representation as well as on defined properties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Event
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".

Function
Function receives argument values, calculates, and returns result of the given expression.

Constraint
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.

Tabular data
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function.
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".
Function law for the reaction, etc.). Based on this visual representation as well as on defined properties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Event
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".

Function
Function receives argument values, calculates, and returns result of the given expression.

Constraint
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.

Tabular data
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function.
Function receives argument values, calculates, and returns result of the given expression.
Constraint law for the reaction, etc.). Based on this visual representation as well as on defined properties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Event
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".

Function
Function receives argument values, calculates, and returns result of the given expression.

Constraint
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.

Tabular data
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function.
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.
Tabular data erties of diagram elements, BioUML automatically generates a program code that is employed to simulate the model dynamics. The current BioUML version generates highly optimized Java code and uses its own state-of-the-art simulation engines.

SBGN Process Diagrams Extension
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 the 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 devised an extension (Table 1) for the Process Diagram type from the SBGN standard [24] to provide the possibility of graphical representation of mathematical elements used in SBML format: equations, events, functions, and constraints [22]. We also added glyphs to represent tabular data that are used for defining the dynamics of the mathematical variables of the model. Tabular data may be translated into spline curves or a constant piecewise function. Furthermore, tabular data may be used, for example, for defining experimental conditions-training regimen for physical exercises.
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports (see below).

Event
Events describing instant changes made to model variables when a specified condition is satisfied, i.e., when trigger expression changes its value from "false" to "true".

Function
Function receives argument values, calculates, and returns result of the given expression.

Constraint
Constraint is checked during the simulation. If it is violated, simulation is either stopped or an error message is shown depending on solver settings.

Tabular data
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function.
Tabular element is used to calculate model variables based on tabular data. In this example, the time column is used for the time variable; x_value and y_value columns are used for x and y, respectively. The tabular element is either translated to a spline approximating tabular data or a piecewise constant function. Interface ports

Interface ports
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports: contact ports (gray), input ports (green), and output ports (red).

Modular Diagrams
Modularity could be considered a principle of biological organization [26,27]. Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last decades [28,29].
In the framework of a modular approach, the investigated system is viewed as a set of interconnected subsystems (modules). Each module is a mathematical model that can be considered and simulated independently. Integration of these modules constitutes a Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports: contact ports (gray), input ports (green), and output ports (red).

Modular Diagrams
Modularity could be considered a principle of biological organization [26,27]. Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last decades [28,29].
In the framework of a modular approach, the investigated system is viewed as a set of interconnected subsystems (modules). Each module is a mathematical model that can be considered and simulated independently. Integration of these modules constitutes a more complex model of the whole system. Modules may leverage different mathematical formalisms and scales. They can be created, validated, and improved independently and may be viewed as replaceable parts.
For visual modelling of modular models, we developed a special diagram type that allows us to specify connections between modules. For this purpose, each module specifies variables that can be used to connect it with other modules. This subset of variables is called the module interface and is represented by ports ( Table 2).
Ports can be of three types: • Input-mathematical variable associated with input ports that is calculated outside of the module and used in the module.

•
Output-mathematical variable associated with contact ports can be modified inside the module as well as outside (e.g., using differential equations).

•
Contact-mathematical variable associated with output ports calculated inside the module and may be used in other modules. In other words, it is a shared variable that can be simultaneously changed by several modules.
Besides this, modular diagrams can include all mathematical elements and tabular data suggested in Table 1. Table 2. Glyphs and arcs for modular diagrams.

Element Name Glyph Description
Module Interface ports Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports: contact ports (gray), input ports (green), and output ports (red).

Modular Diagrams
Modularity could be considered a principle of biological organization [26,27]. Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last decades [28,29].
In the framework of a modular approach, the investigated system is viewed as a set of interconnected subsystems (modules). Each module is a mathematical model that can be considered and simulated independently. Integration of these modules constitutes a more complex model of the whole system. Modules may leverage different mathematical formalisms and scales. They can be created, validated, and improved independently and may be viewed as replaceable parts.
For visual modelling of modular models, we developed a special diagram type that allows us to specify connections between modules. For this purpose, each module specifies variables that can be used to connect it with other modules. This subset of variables is called the module interface and is represented by ports ( Table 2).
Ports can be of three types: • Input-mathematical variable associated with input ports that is calculated outside of the module and used in the module.

•
Output-mathematical variable associated with contact ports can be modified inside the module as well as outside (e.g., using differential equations).

•
Contact-mathematical variable associated with output ports calculated inside the module and may be used in other modules. In other words, it is a shared variable that can be simultaneously changed by several modules.
Besides this, modular diagrams can include all mathematical elements and tabular data suggested in Table 1. Table 2. Glyphs and arcs for modular diagrams.

Module
Module encapsulates the mathematical model (submodel) of a particular subsystem forming the hierarchic structure of the model.

Port
Ports represent the interface of the module through which it can be connected with other modules or with the hierarchical model itself. The color of the port defines its type: output (red), input (green), contact (gray).

Connection
Connections can be established between ports to aggregate modules with each other. Directed connections are established between output and input ports, with undirected connections between contact ports.
Module encapsulates the mathematical model (submodel) of a particular subsystem forming the hierarchic structure of the model.

Interface ports
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports: contact ports (gray), input ports (green), and output ports (red).

Modular Diagrams
Modularity could be considered a principle of biological organization [26,27]. Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last decades [28,29].
In the framework of a modular approach, the investigated system is viewed as a set of interconnected subsystems (modules). Each module is a mathematical model that can be considered and simulated independently. Integration of these modules constitutes a more complex model of the whole system. Modules may leverage different mathematical formalisms and scales. They can be created, validated, and improved independently and may be viewed as replaceable parts.
For visual modelling of modular models, we developed a special diagram type that allows us to specify connections between modules. For this purpose, each module specifies variables that can be used to connect it with other modules. This subset of variables is called the module interface and is represented by ports ( Table 2).
Ports can be of three types: • Input-mathematical variable associated with input ports that is calculated outside of the module and used in the module.

•
Output-mathematical variable associated with contact ports can be modified inside the module as well as outside (e.g., using differential equations).

•
Contact-mathematical variable associated with output ports calculated inside the module and may be used in other modules. In other words, it is a shared variable that can be simultaneously changed by several modules.
Besides this, modular diagrams can include all mathematical elements and tabular data suggested in Table 1. Table 2. Glyphs and arcs for modular diagrams.

Module
Module encapsulates the mathematical model (submodel) of a particular subsystem forming the hierarchic structure of the model.

Port
Ports represent the interface of the module through which it can be connected with other modules or with the hierarchical model itself. The color of the port defines its type: output (red), input (green), contact (gray).

Connection
Connections can be established between ports to aggregate modules with each other. Directed connections are established between output and input ports, with undirected connections between contact ports.
Ports represent the interface of the module through which it can be connected with other modules or with the hierarchical model itself.
The color of the port defines its type: output (red), input (green), contact (gray).

Interface ports
Although SBGN notation already has tag elements that denote the module interface (ports in SBML terminology), in our diagrams we have three different types of ports: contact ports (gray), input ports (green), and output ports (red).

Modular Diagrams
Modularity could be considered a principle of biological organization [26,27]. Therefore, a modular approach to the modeling of complex biochemical systems has been actively developing in the last decades [28,29].
In the framework of a modular approach, the investigated system is viewed as a set of interconnected subsystems (modules). Each module is a mathematical model that can be considered and simulated independently. Integration of these modules constitutes a more complex model of the whole system. Modules may leverage different mathematical formalisms and scales. They can be created, validated, and improved independently and may be viewed as replaceable parts.
For visual modelling of modular models, we developed a special diagram type that allows us to specify connections between modules. For this purpose, each module specifies variables that can be used to connect it with other modules. This subset of variables is called the module interface and is represented by ports ( Table 2).
Ports can be of three types: • Input-mathematical variable associated with input ports that is calculated outside of the module and used in the module.

•
Output-mathematical variable associated with contact ports can be modified inside the module as well as outside (e.g., using differential equations).

•
Contact-mathematical variable associated with output ports calculated inside the module and may be used in other modules. In other words, it is a shared variable that can be simultaneously changed by several modules.
Besides this, modular diagrams can include all mathematical elements and tabular data suggested in Table 1. Table 2. Glyphs and arcs for modular diagrams.

Module
Module encapsulates the mathematical model (submodel) of a particular subsystem forming the hierarchic structure of the model.

Port
Ports represent the interface of the module through which it can be connected with other modules or with the hierarchical model itself. The color of the port defines its type: output (red), input (green), contact (gray).

Connection
Connections can be established between ports to aggregate modules with each other. Directed connections are established between output and input ports, with undirected connections between contact ports.

Visual Modular Modelling
Module ports are used on two levels ( Figure 1): • when creating model that will be used as a part of another model (i.e., module), a modeler specifies the module inputs, outputs, and contacts and links them to corresponding module entities or variables; • on a modular diagram, a modeler links several modules together using previously defined ports.
Let us consider a simple example demonstrating this approach ( Figure 1). First, we will develop a simple module M1 that consists of one biochemical reaction where two molecules A and B form the complex A:B. We are suggesting that the concentration of entity A can be changed in other reactions due to participation in other modules. To spec-Buses are auxiliary elements that can be used as intermediate points in connections. They decrease the number of intersections between connections and make the structure of the model more clear.

Visual Modular Modelling
Module ports are used on two levels ( Figure 1): • when creating model that will be used as a part of another model (i.e., module), a modeler specifies the module inputs, outputs, and contacts and links them to corresponding module entities or variables; • on a modular diagram, a modeler links several modules together using previously defined ports.
Let us consider a simple example demonstrating this approach ( Figure 1). First, we will develop a simple module M1 that consists of one biochemical reaction where two molecules A and B form the complex A:B. We are suggesting that the concentration of entity A can be changed in other reactions due to participation in other modules. To specify this, we will create port A of type "contact" (grey pentagon). The concentration of the A:B complex is solely defined in module M1, and we will create a port A:B of type "output" (red pentagon) that is represented as input in module M2.
Then, we will create module M2 that will also include one reaction where complex A:B catalyzes the phosphorylation of protein C (C{p}). Herein, we will define port A:B as input (green pentagon) for this module and port C{p} as output.
Module M3 also comprises one bimolecular reaction where C{p} catalyzes the transformation of A into A_1. X and Z are chemical substances that are the reactant and product, respectively. Similarly, we will specify C{p} as input, and A port will be a contact while the A concentration is also changed in the reaction from module M1.
Finally, let us form these three modules into a modular model (Figure 1b,c). We will connect corresponding ports to each other ( Figure 1b).
• inputs and outputs: A:B for M1 and M2, C{p} for M2 and M3; this is a directed connection so it is depicted by an arrow; • contacts-A for M1 and M3; this is an undirected connection while concentration A is changed simultaneously by two reactions from these modules and so it is depicted by the line.  . M1-M3 designate corresponding modules; two green rectangles A and B in M1 correspond to molecules A and B which form the complex A:B, while grey and red pentagons in M1 designate contact port for A and output port for A:B, respectively; two green rectangles C and C_p in M2 correspond to protein C and phosphorylated form of the protein, while green and red pentagons in M2 designate input port for A:B and output port for C_p, respectively; two green rectangles A and A_1 in M3 correspond to molecules A and A_1, while green and grey pentagons in M3 designate input port for C_p and contact port for A, respectively; two purple circles X and Z in M3 mean the additional substrate and product of the bimolecular reaction, correspondingly that is catalyzed by phosphorylated form of the protein C (green rectangle), (b,c)-modular diagram in two equivalent variants: without or with buses (white circles). Buses serve for cosmetic purposes only.
Numerical calculations for modular models may be performed in two ways: • Flattening-a modular model may be transformed into a non-modular model by aggregating all elements of all modules with automatic resolving of established connections between variables [18]. which form the complex A:B, while grey and red pentagons in M1 designate contact port for A and output port for A:B, respectively; two green rectangles C and C_p in M2 correspond to protein C and phosphorylated form of the protein, while green and red pentagons in M2 designate input port for A:B and output port for C_p, respectively; two green rectangles A and A_1 in M3 correspond to molecules A and A_1, while green and grey pentagons in M3 designate input port for C_p and contact port for A, respectively; two purple circles X and Z in M3 mean the additional substrate and product of the bimolecular reaction, correspondingly that is catalyzed by phosphorylated form of the protein C (green rectangle), (b,c)-modular diagram in two equivalent variants: without or with buses (white circles). Buses serve for cosmetic purposes only.
More complex modular diagrams may contain a large number of connections that form many intersections. To skip this intersection, we introduce the concept of a bus (white circle in Figure 1c): a port is connected to the named bus, and a diagram may contain several clones of such a bus. Figure 1c demonstrates how the connection of two A ports can be replaced by a connection with two clones of bus A.
Numerical calculations for modular models may be performed in two ways: • Flattening-a modular model may be transformed into a non-modular model by aggregating all elements of all modules with automatic resolving of established connections between variables [18]. • Agent-based simulation. Each module is simulated independently with its own simulator and formalism. The implemented scheduler coordinates the interactions by sending and receiving the numerical values of connected variables [19].
To simulate the presented integrated modular model, we employ a flattening approach while all modules use the same mathematical formalism and contain sets of ordinary differential equations (ODEs) and discrete events (i.e., hybrid models). The BioUML platform automatically transforms the modular model into a "flat" hybrid model with the same formalism by aggregating all equations and events from all modules and resolving connections. For more details, see [18,19].

Antimony-Extension and Synchronization with Visual Depiction
Antimony provides a convenient human-readable text format that supports most of the SBML features. Herein, we suggest an extension for the Antimony format to specify reaction components to which type of SBGN Process Diagram they correspond, as well as some other properties.
The suggested format is as follows: @entity_id.property = value The proposed extension is quite similar to the idea of annotations in SBML format where SBML-enabled software can store any auxiliary information. Similar to Java annotations, we suggest sign @ for this purpose.
Currently, the BioUML platform supports the following properties in annotations: • sbgnType-defines the SBGN entity type (unspecified, macromolecule, nucleic acid feature, perturbing agent, simple chemical or complex). All those entities correspond to mathematical variables in the model. • sbgnViewTitle-defines additional properties of an SBGN entity such as whether it is multimer if it has units of information or state variables. If an entity is a complex, it also defines the inner elements of the said complex. We used Transpath conventions to denote entities and complexes in text formats. Here are some examples: Complex comprising two entities A and B is denoted as "A:B". Entity A with state variable p (phosphorylated) is denoted as "A{p}". Multimer entity A is denoted as "(A)3". A more advanced example: "(A{p})3:B{r}{p}:C".
Depending on the context/tasks, it can be more suitable to present a model of a biological system either as a graph using the extended version of the SBGN Process diagram type or as a program code using Antimony language.
Antimony+ and PD+ are seamlessly integrated in the frame of the BioUML platform. Due to this integration, a user can simultaneously view and edit textual and graphical representations of a biological system model. Figure 2 demonstrates how the chemical reaction is represented using SBGN Process Diagram Type (2a) and extended Antimony format (2b).
logical system either as a graph using the extended version of the SBGN Process diagram type or as a program code using Antimony language. Antimony+ and PD+ are seamlessly integrated in the frame of the BioUML platform. Due to this integration, a user can simultaneously view and edit textual and graphical representations of a biological system model. Figure 2 demonstrates how the chemical reaction is represented using SBGN Process Diagram Type (2a) and extended Antimony format (2b). It is noteworthy that if a user edits textual model representation then graphical representation is updated synchronously by the BioUML platform and vice versa. Similarly, if a user selects some object on a diagram, then corresponding text items are highlighted in the text document and vice versa. It is noteworthy that if a user edits textual model representation then graphical representation is updated synchronously by the BioUML platform and vice versa. Similarly, if a user selects some object on a diagram, then corresponding text items are highlighted in the text document and vice versa.

Model Simulations
Numerical solutions of the model represented by a system of ordinary differential equations have been obtained on the basis of the VODE method [30] using a JVode simulation engine implemented in the BioUML tool [17]. Each submodule of the modular model can be represented as an independent SBML file [22], while the integrated modular model can be exported as a COMBINE archive [31] to use the model and reproduce simulations, resulting in alternative software supporting current standards of the systems biology.

Jupyter Notebook
BioUML is integrated with Jupyter (https://jupyter.org, accessed on 8 June 2021) for interactive data and model analysis as well as an essential and user-friendly tool for the reproducibility of the simulation results ( Figure 3). The notebook that provides reproducibility of results presented in the article can be started on the BioUML server as well as using Binder technology.

The Model Overview
The general structure of the developed model linking metabolism, Ca 2+ -dependent signaling transduction, and regulation of gene expression in human skeletal muscle is demonstrated in Figure 4.
The model has a hierarchical structure. At the top level, the model describes the physiology of capillary blood flow through muscles during exercise to provide oxygen and substrate delivery and metabolite removal from the skeletal muscle. It is worth noting that many physiological details are lumped in the current version of the model (e.g., cardiorespiratory system elements), and the dynamic change of the capillary blood flow elicited by the exercise is considered a linear function of the exercise intensity [32]. In the muscle model, we consider Type I and Type II fibers. Their modules have the same structure but differ in some parameter values.

Jupyter Notebook
BioUML is integrated with Jupyter (https://jupyter.org, accessed on 24 September 2021) for interactive data and model analysis as well as an essential and user-friendly tool for the reproducibility of the simulation results ( Figure 3). The notebook that provides reproducibility of results presented in the article can be started on the BioUML server as well as using Binder technology.

The Model Overview
The general structure of the developed model linking metabolism, Ca 2+ -dependent signaling transduction, and regulation of gene expression in human skeletal muscle is demonstrated in Figure 4.
The model has a hierarchical structure. At the top level, the model describes the physiology of capillary blood flow through muscles during exercise to provide oxygen and substrate delivery and metabolite removal from the skeletal muscle. It is worth noting that many physiological details are lumped in the current version of the model (e.g., cardiorespiratory system elements), and the dynamic change of the capillary blood flow elicited by the exercise is considered a linear function of the exercise intensity [32]. In the muscle model, we consider Type I and Type II fibers. Their modules have the same structure but differ in some parameter values. On the next level (cellular level), we consider biological processes that occur in human skeletal muscle cells. These processes can also be divided into three sublevels that are described by the corresponding modules of the model:

•
Metabolic-the main metabolic processes that occur in the skeletal muscle during physical exercises: glycolysis, glycogenolysis, tricarboxylic acid cycle, β-oxidation, and oxidative phosphorylation. This part of the model is based on a detailed mathematical model of muscle metabolism developed by Li and coauthors [32]. We have redesigned this model according to the methodology described above and changed some model parameters to reproduce more experimental data (see below). • Signalling-the main signal transduction pathways that are activated by physical exercises are related to Ca 2+ -dependent signaling [33] and AMPK activation [34]. For each of them we developed a special module. • Gene expression regulation-changes in gene expression were divided into early (up to 1-3 h after exercises) and late (3-6 h after exercises) responses. We selected the most well-known genes for each group-NR4A2 and NR4A3 for the first group and PPRGC1A for the second as described in the "Gene expression level" section. The corresponding modules that describe the expression of these genes have also been developed.
Oxygen delivery and metabolite transport between cellular compartments (mitochondria and cytoplasm) as well as between muscle cells and capillary blood are very important parts of the model. Therefore, we developed special modules considering all these transport processes.
An activation mechanism that enhances energy metabolism via transport and reaction fluxes due to physical exercise was harnessed as the stress function depending on the general work rate parameter [32]: where α i is the activation coefficient, W is the work rate value, τ i indicates the time constant of changes in the metabolic reaction rates in response to exercise, and t start is the simulation time when an exercise is started. The work rate parameter defines the power of the physical exercise and varies depending on the mode of the exercise. In our model, the muscle volume was 5 L, which corresponds to that involved in exercise using a cycling ergometer [14,32].   All details on the general description of each module, representation of the module diagram using extended SBGN Process Description notation, corresponding Antimony code for the module as well as reaction rates and ordinary differential equations to describe the species concentrations, the module parameters, their values, and references to the literature from which they were extracted are presented in the Supplementary Materials. Overall, the model includes: • 25 modules; • 238 species; • 185 reactions; • 171 ordinary differential equations; • 647 parameters.

Physiological (or Organism) Level
On this level ( Figure 5), we model capillary blood flow for oxygen and substrate delivery to the muscle and for removal of metabolites produced by the muscle including: CO 2carbon dioxide; O 2 -oxygen; Lac-lactate; Ala-alanine; Pyr-pyruvate; H-hydrogen; Glr-glycerol; Glc-glucose; FFA-free fatty acid. The model comprises five modules: "Capillary_Blood_Flow" to describe flow processes in the blood; "Cytosol_Capillary Transport R" and "Cytosol_Capillary Transport W" modules representing transport reactions between capillary blood and the muscle tissue, "Fiber R" and "Fiber W" modules where metabolic, signaling, and gene expression regulation processes are considered. All details for each submodule and zoomed-in subfigures are in the Supplementary Material.

Transport Level
The inter-compartmental metabolite transport is described as passive or facilitated (carrier mediated) fluxes according to the original paper [32]. By analogy with metabolic rates, all transport flux equations are multiplied by a linear function to consider the exercise effect on transport processes. The basic transport flux equation for passive (superscript p) diffusion of species i between the blood and cytosol is: Figure 5. The modular representation of the integrated model at the physiological level. The model comprises five modules: "Capillary_Blood_Flow" to describe flow processes in the blood; "Cy-tosol_Capillary Transport R" and "Cytosol_Capillary Transport W" modules representing transport reactions between capillary blood and the muscle tissue, "Fiber R" and "Fiber W" modules where metabolic, signaling, and gene expression regulation processes are considered. All details for each submodule and zoomed-in subfigures are in the Supplementary Material.
It should be noted that those species are present in six different modules and have different subscripts. We explain this using the example of CO 2 . In the module "Capil-lary_Blood_Flow", CO 2c is the concentration of CO 2 in the capillary blood. It is connected via connections and buses with CO 2b in modules "Cytosol_Capillary Transport R" and "Cytosol_Capillary Transport W", where CO 2b is also the CO 2 concentration in the capillary blood. In those modules, transport of CO 2 from blood to the muscle tissue is presented, where CO 2c is the concentration in the tissue. Finally, CO 2c is connected with the CO 2 variables in the "Fiber R" and "Fiber W" modules (corresponding to Type I and Type II fibers, respectively), where the metabolism of CO 2 in tissues is considered. Bus elements (white circles in Figure 5) are used to prevent too many intersections between connections.
Skeletal muscle volume (V mus ) is represented by the sum of the effective volumes of blood (V bl ): V mus = V tis + V bl -the skeletal muscle volume (5 kg w.w.). In skeletal muscle, recruitment of muscle fibers due to the transition from rest to exercise enhances active muscle volume and blood flow. According to the original metabolic model [32], we also assume that these physiological variables dynamically change as linear functions of the work rate parameter or power of the physical exercise: where Q is the blood flow, Q 0 = 0.9 L/min is the muscle blood flow at rest for two legs; V mus 0 is the skeletal muscle volume at rest (5 kg w.w.); α i is the activation coefficient; while τ V = τ Q = 0.4 min is the time constant of the muscle volume and blood flow changes, respectively, in response to exercise; and t start is the simulation time when an exercise is started [32]. The muscle consists of two compartments (modules) that correspond to Type I and Type II fibers. They have the same structure but differ in some parameters (see details in corresponding tables of the Supplementary Materials).

Transport Level
The inter-compartmental metabolite transport is described as passive or facilitated (carrier mediated) fluxes according to the original paper [32]. By analogy with metabolic rates, all transport flux equations are multiplied by a linear function to consider the exercise effect on transport processes. The basic transport flux equation for passive (superscript p) diffusion of species i between the blood and cytosol is: where λ bl − cyt,type, i is the permeability-surface area coefficient, C bl, i and C cyt,type,i are concentrations of species i in the blood and cytosol, respectively; i ∈ (CO 2 , O 2 , Ala, Glr) and type ∈ (type I f iber, type I I f iber), while for facilitated (superscript f ) transport: where R max transport cyt − mit, type,i is the maximal flux rate for facilitated transport, C bl, i and C cyt,type,i are concentrations of species i in the blood and cytosol, respectively; i ∈ (Glc, Pyr, Lac, FFA, H + ) and type ∈ (type I f iber, type I I f iber).
The basic transport flux equation for passive (superscript p) diffusion of species i between the cytosol and mitochondria is: where λ cyt − mit,type, i is the permeability-surface area coefficient, C cyt,type, i and C mit,type,i are concentrations of species i in cytosol and mitochondria, respectively; i ∈ (CO 2 , O 2 ) and type ∈ (type I f iber, type I I f iber), while for facilitated (superscript f ) transport: where R max transport cyt − mit, type,i is the maximal flux rate for facilitated transport, C cyt,type, i and C mit,type,i are concentrations of species i in the cytosol and mitochondria, respectively; i ∈ (H + , Pyr, FAC, CoA, P i ) and type ∈ (type I f iber, type I I f iber).

Cellular (Metabolic) Level
The diagram of the modular model describing the metabolism of human skeletal muscle is presented in Figure 6. The cytosol includes metabolic reactions of the glycolysis, glycogenolysis, and lipid metabolism, while the tricarboxylic acid (TCA) cycle, ß-oxidation, and oxidative phosphorylation reactions occur 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 coauthors [32]. According to the model, a dynamic mass balance of metabolites (i) is based on the structural and functional organization and is expressed by equations: where V cyt , V mit indicate the volume of the corresponding module or compartment in kg wet weight (kg w.w.), type ∈ (type I f iber, type I I f iber). V cyt,R = 0.88 × V R and V cyt,W = 0.92 × V W are volumes of the cytosol for type I and II fibers, respectively, while V mit,R = 0.12 × V R and V mit,W = 0.08 × V W are the volumes of mitochondria for type I and II fibers, respectively, where V R = V w = 0.5 × V tis = 2 kg w.w., V R -the type I fiber volume, V W -the type II fiber volume, and V tis -the volume of muscle cells in the tissue. C cyt, type,i , C mit, type,i is the concentration of metabolite i in a certain compartment of the corresponding fiber type (mmol/kg w.w.); R x,type,i , x ∈ {cyt, mit} is the net metabolic reaction rate in a certain compartment of the corresponding fiber type (mmol/min/kg w.w.); T k cyt − mit, type,i , T k bl − cyt, type,i are the respective transport fluxes between the cytosol and mitochondria compartments and cytosol and blood compartments (mmol/kg w.w.), where superscript k indicates f (facilitated) or p (passive) transports.
In order to describe a dynamic mass balance of metabolites (i) in the blood compartment, the next equation is used: where V bl is the total effective volume of the capillary blood and interstitial fluid compartments. V bl = 0.2 × V mus , V tis = 0.8 × V mus , where V mus = V tis + V bl -the skeletal muscle volume (5 kg w.w.); C art,i C bl,i is the concentration of metabolite i in the respective arterial and capillary blood compartments (mmol/kg w.w.). It is worth noting that such modules as capillary blood and interstitial fluid are assumed to be in equilibrium with each other, which allows us to consider species concentrations in these compartments as equal and combine them into the blood compartment. The comprehensive description including visual representation, corresponding Antimony code, and mathematical equations for reaction rates and the dynamic mass balance in each module of the integrated model as well as values and units of the used kinetic parameters is presented in the Supplementary material.
In comparison with the original model of Li and coauthors [32], we introduced a following changes:

1.
Values of activation coefficients associated with ATPase [35][36][37][38] and pyruvate dehydrogenase reaction fluxes for type I and type II fibers [39][40][41] as well as the time constant of the ATPase flux rate coefficient in response to exercise were modified (See Data availability and Supplementary material) according to recently published data and estimations [42,43].

2.
Despite overall net glycogen breakdowns during muscle contraction, exercise increases the activity of glycogen synthase (GS) [44][45][46][47] and ATP consumption related with the reaction. Therefore, GS reaction fluxes were modified according to [44,46,48]. The rates of muscle glycogen synthesis during exercise assumed to be equal in type I and type II fibres were estimated from average post-exercise glycogen synthesis data [49].

3.
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 [50][51][52].

Signaling Level
The mean concentration of Ca 2+ ions in the myoplasm increases in proportion to the intensity of exercise. Ca 2+ binds to calmodulin, thereby activating CaMKs and phosphatase calcineurin [33]. CaMKII is the most abundant isoform in human skeletal muscle, whereas CaMKI and CaMKIV are not expressed at detectable levels [53]. An increase in CaMKII activity results in CREB1 Ser133 phosphorylation (and likely some other CREB-like factors), leading to the activation of the transcription factor [54,55]. Calcineurin can dephosphorylate (and activate) CRTCs at Ser171 (CREB-regulated transcription coactivators), playing a key role in regulating the transcriptional activity of CREB1 [56]. Another target of calmodulin is Ca 2+ /calmodulin-dependent protein kinase kinase 2 (CAMKK2) that phosphorylates AMPK Thr172, thereby activating the kinase [57]. In turn, activated AMPK can phosphorylate CREB1 Ser133 [58]. Collectively, these findings drove us to include in our model the Ca 2+ -dependent regulation of calmodulin, CREB1 (via CaMKII), CRTC (via calcineurin), and AMPK (via CaMKK2) (Figure 7). The amount of these proteins in human skeletal muscle was estimated using published proteomics and transcriptomics data [12,59] (see Supplementary data in [60]).
There are three different heterotrimeric complexes of AMPK in the human skeletal muscles: α2β2γ1, α2β2γ3, and α1β2γ1 [61]. Distinct kinetic properties (an intrinsic enzyme activity, binding affinities of AMP, ADP, and ATP to the specific isoform, sensitivity to deand phosphorylation of AMPK heterotrimers) [62,63] and their subcellular localization [64] cause a differential regulation of the AMPK heterotrimers in vivo. The α2β2γ3 complex is phosphorylated and activated during moderate-to high-intensity aerobic exercise, while the activity associated with the other two AMPK heterotrimers is almost unchanged [65]. However, the baseline activity of the α2β2γ3 complex is significantly lower than others. The general AMPK activity at baseline and during/after exercise is a sum of isoform activities; hence, we considered all isoforms separately (in the corresponding module) to quantitatively fit experimental data obtained at baseline and after an exercise [65,66]. AMPK is regulated in various ways: an up-stream kinase LKB1 can phosphorylate AMPK at Thr172 [67,68]. On the other hand, an exercise-induced decrease in intramuscular ATP increases its activity, while an increase in AMP activates it [69,70]. Hence, in our model, AMPK is regulated via AMP, ATP, and LKB1, as well as CaMKK2 (as mentioned above) (Figure 7). isoform activities; hence, we considered all isoforms separately (in the corresponding module) to quantitatively fit experimental data obtained at baseline and after an exercise [65,66]. AMPK is regulated in various ways: an up-stream kinase LKB1 can phosphorylate AMPK at Thr172 [67,68]. On the other hand, an exercise-induced decrease in intramuscular ATP increases its activity, while an increase in AMP activates it [69,70]. Hence, in our model, AMPK is regulated via AMP, ATP, and LKB1, as well as CaMKK2 (as mentioned above) (Figure 7).   It was demonstrated that the localizations of AMPK and CaMKII kinases have a pronounced effect on their activities [34,53,[71][72][73], implying the necessity to consider the impact in the model. However, an extended analysis of the published data on this issue demonstrates some contradictions in the data and the lack of quantitative data on this issue. For instance, the vast majority of CaMKII (~80%) expressed in human skeletal muscle is localized to the soluble cytosolic fraction. However, most of the major estimations and measurements on the functional properties and substrates have been obtained for membrane-associated CaMKII [53]. Moreover, the mobile fraction of the kinases or their substrates has a limited diffusion rate in the tightly packed myocyte structure and is dependent on the molecular weight that can affect the kinetics of their interaction. Such diffusion rate data have not been found.

Gene Expression Level
An aerobic exercise induces the expression of several hundreds of genes regulating many cell functions: energy metabolism, transport of various substances, angiogenesis, mitochondrial biogenesis, etc. Regulation of the transcriptomic response to acute exercise includes dozens of transcription regulators [12] and seems to be extremely complex. Therefore, to consider the response at a gene expression level, we selected exercise-induced genes based on the next criterion comprising two points: (1) a functional role of this gene in the regulation of skeletal muscle metabolism is known; (2) its expression in human skeletal muscle markedly changes in response to an exercise and has one of the expression patternsearly or delayed response since gene expression in early and late stages of the recovery after the termination of the exercise can be regulated in different ways [13]. According to the criterion, the PPARGC1A gene, known as the major regulator of exercise-induced phenotypic adaptation and substrate utilization [74], was chosen as the gene with delayed response, while NR4A2 and NR4A3 genes were chosen as early response genes [75]. NR4A nuclear receptors induce DNA demethylation in skeletal muscle [76], regulate genes involved in glycogenolysis, fatty acid oxidation, and pyruvate use and apparently play a role in the regulation of the skeletal muscle fiber phenotype [77,78]. Significantly, all members of the NR4A nuclear receptor subfamily (NR4A1, NR4A2, NR4A3) are the three most highly induced genes in response to acute endurance exercise [79,80]. We selected both genes from one family since they have different temporal patterns of mRNA expression that are likely associated with different methylation profiles of their promoters [81,82].
Expression of NR4A2 and NR4A3 mRNA rapidly increases during the first hour after an aerobic exercise (early response genes) [12] due to activation of Ca 2+ \calcineurindependent signaling [75]. We included in our model the Ca 2+ -dependent regulation (Ca 2+ \calcineurin-CaMKII-CREB1) of NR4As genes using data of contractile activityspecific mRNA responses of these genes [12]. Expression of PPARGC1A mRNA rises 3 to 4 h after an exercise (gene with delayed response) [12]. The transcription regulation of PPARGC1A via the canonical (proximal) and inducible (distal) promoters is very complicated and includes Ca 2+ -and AMPK-dependent signaling, as well as CREB1 and its co-activator CRTC [10,83]. The phosphorylation level of many signaling kinases drops to basal levels within the first hour after an aerobic exercise. Moreover, in a genomewide study on various human tissues, it was shown that the phosphorylation level of CREB Ser133 does not always correlate with its transcriptional activity [81]. Therefore, we suggested that the expression of genes with delayed response (including PPARGC1A) is regulated by increasing the expression of one of the early response genes encoding transcription factors leading to a rapid increase in the corresponding protein [60]. A detailed description of our results on the identification of transcription factors as potential candidates for the role of X factor is presented below in the section Results and Discussion. Analysis of contractile activity-specific transcriptomic data [12] 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 intersect with each other both in the alternative and in the canonical promoters of the PPARGC1A gene [60], i.e., these TFs can act as 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 homoand heterodimers [84,85]. 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 existing TFs (e.g., CREB1) and their co-activators (e.g., CRTC), while delayed response genes are regulated via an increase in the expression of early response genes encoding transcription factors (transcription factor X in our model, Supplementary Material, Module "Gene expression regulation").

Simulation of Metabolic Changes Induced by Incremental and Interval Exercises
To validate the metabolic part of the model, we investigated the dynamic behaviour of the system in response to diverse acute aerobic exercises and compared them with published experimental data. It is worth noting that qualitative validation of the model was conducted on the basis of the comparison of the simulation and experimental data for three indicators: time period to achieve the maximal level of the species concentrations (e.g., PCr, ATP, glycogen) at the corresponding value of the exercise intensity and time to reach the steady-state value in recovery as well as the multiplicity of concentration changes (fold changes). We used the last indicator due to quantitative differences in measured concentrations for the same species by different experimental approaches (e.g., biochemical and 31 P MRS measurements). Initially, we quantitatively estimated the biochemical responses of the key metabolic variables (ATP, ADP, PCr, lactate concentrations, and pH in muscle fibers type I and II) in the incremental ramp exercise to exhaustion, which is a commonly used approach to evaluate aerobic performance. Increasing the power during the ramp exercise affects various physiological variables such as the number/volume of recruited muscle fibre type I and II, blood flow as well as the transport and metabolic fluxes in both fibre types (Figure 8 and see data availability). In our simulation, muscle fibres type I start to be recruited after the beginning of exercise, while fibre type II if only recruited at a power higher than 24% of VO2 max (6 min after the ramp exercise onset, Figure 8A). Recruiting all muscle fibres during the test leads to exhaustion and termination of the exercise [35,86,87]; the peak power at exhaustion in our simulation was 250 W, which corresponds to the value for an untrained male performing the ramp exercise until exhaustion using a cycling ergometer. The model simulations correspond reasonably well to experimental measurements [88][89][90] obtained in studies with the incremental exercise (Supplementary Figure S1). It is worth noting that the current version of the model does not take into account the effect of muscle fatigue during the incremental ramp exercise observed in exercised muscle in vivo (see below). This fact may partially explain the lack of exponential changes in muscle lactate concentration and pH during the last part of the incremental exercise.  For additional validation of the metabolic part of the model, we simulated responses to various interval exercises (Figures 9 and 10). Figure 9 shows that the model qualitatively reproduces the dynamics of PCr concentration during the interval exercises with different patterns (duration of an exercise bout 16 s to 64 s and recovery period 32 s to 128 s) and with peak power comparable with maximal aerobic power obtained in the incremental ramp test (250 W) [91]. For additional validation of the metabolic part of the model, we simulated responses to various interval exercises (Figures 9 and 10). Figure 9 shows that the model qualitatively reproduces the dynamics of PCr concentration during the interval exercises with different patterns (duration of an exercise bout 16 s to 64 s and recovery period 32 s to 128 s) and with peak power comparable with maximal aerobic power obtained in the incremental ramp test (250 W) [91].
For additional validation of the metabolic part of the model, we simulated responses to various interval exercises (Figures 9 and 10). Figure 9 shows that the model qualitatively reproduces the dynamics of PCr concentration during the interval exercises with different patterns (duration of an exercise bout 16 s to 64 s and recovery period 32 s to 128 s) and with peak power comparable with maximal aerobic power obtained in the incremental ramp test (250 W) [91]. High-intensity interval exercise has been shown to rapidly decrease the PCr level followed by slow recovery of the PCr concentration during the last part of the exercise [92]. There are no data on the exercise power in the study; hence, we used the constant value (500 W) for each bout ( Figure 10A). The power was markedly higher than the peak power in the incremental cycling test because the duration of each exercise bout is short; the energy supply of such short exercise bouts is related mainly to PCr reactions as well as glycolysis. Our model precisely simulated the rapid decline in PCr, but showed no slow recovery of the PCr during the last part of the high-intensity interval exercise ( Figure 10B). We suggested that this discrepancy may be related to the lack of the fatigue-induced decline in exercise power. We tried to roughly simulate the fatigue-induced decline in exercise power by the decline in power generated by muscle fibers type II ( Figure 10C). As a result, the model much better reproduced the experimental dynamics of PCr than simulations with constant maximal power in each bout ( Figure 10D,E, Supplementary Figure S2). However, the PCr dynamics during the recovery process indicated that the model still requires further modifications and numerical study. We assume that the potential point for the update is related to the pH changes during the recovery.

Simulation of Signaling and Gene Expression Changes Induced by Low-and Moderate Intensity Continuous Exercises
At the next step of the model validation, we predicted the responses of biochemical variables, signaling molecules (AMPK and Ca 2+ -dependent proteins), transcription factor (CREB1), as well as expression of genes with early and delayed responses (NR4A3, NR4A2, PPARGC1A) to low (50% VO2 max ) and moderate intensity (70% VO2 max ) continuous aerobic exercises ( Figure 11). Moderate intensity exercise recruits more muscle fibers type II than low intensity exercise, thereby additionally modulating the exercise-induced metabolic fluxes and molecular response. A comparison of our simulations with experimental data [43,90,[93][94][95] showed that the model well reproduces the metabolic changes in various fiber types and in the whole muscle induced by exercises with various intensity (Supplementary Figure S1).
According to the literature data on the human vastus lateralis muscle [53,65,96], our simulation showed an intensity-dependent increase in the phosphorylation of CAMKII and AMPK α2 and γ3 ( Figure 12A-C,F-H). Importantly, the phosphorylation (as a marker of activity) of AMPK α2 and γ3 consisted of 10% and 30% of the AMPK isoforms containing α2 and γ3, respectively (Figure 12B-C,G-H), which is in line with the experimental data [65]. In contrast to experimental data at a signaling level, we found transcriptomics data concerning intensity-dependent gene expression for 1 h exercise only [10,12]. In our model, exerciseinduced activation of CAMKII and AMPK induced CREB-and CRTC-related expression of early response genes that is in line with the experimental data [12] on exercise-induced expression of early response genes (for example, of NR4A2, NR4A3) in the human vastus lateralis muscle ( Figure 12K). value (500 W) for each bout ( Figure 10A). The power was markedly higher than the peak power in the incremental cycling test because the duration of each exercise bout is short; the energy supply of such short exercise bouts is related mainly to PCr reactions as well as glycolysis. Our model precisely simulated the rapid decline in PCr, but showed no slow recovery of the PCr during the last part of the high-intensity interval exercise ( Figure 10B). We suggested that this discrepancy may be related to the lack of the fatigue-induced decline in exercise power. We tried to roughly simulate the fatigue-induced decline in exercise power by the decline in power generated by muscle fibers type II ( Figure 10C). As a result, the model much better reproduced the experimental dynamics of PCr than simulations with constant maximal power in each bout ( Figure 10D,E, Supplementary Figure  S2). However, the PCr dynamics during the recovery process indicated that the model still requires further modifications and numerical study. We assume that the potential point for the update is related to the pH changes during the recovery.

Simulation of Signaling and Gene Expression Changes Induced by Low-and Moderate Intensity Continuous Exercises
At the next step of the model validation, we predicted the responses of biochemical variables, signaling molecules (AMPK and Ca 2+ -dependent proteins), transcription factor (CREB1), as well as expression of genes with early and delayed responses (NR4A3, NR4A2, PPARGC1A) to low (50% VO2max) and moderate intensity (70% VO2max) continuous aerobic exercises ( Figure 11). Moderate intensity exercise recruits more muscle fibers type II than low intensity exercise, thereby additionally modulating the exercise-induced metabolic fluxes and molecular response. A comparison of our simulations with experimental data [43,90,[93][94][95] showed that the model well reproduces the metabolic changes in various fiber types and in the whole muscle induced by exercises with various intensity (Supplementary Figure S1).  [92]). (A,C) Exercise power and fiber recruitment pattern: total power (A: W peak = 500; C: W peak = 500), W (orange), power generated by type I (red, A-C: W peak = 125) and II (blue, A: W peak = 375; C: W peak = 375 and successive power decline) fibers; (B,D) PCr concentration in type I (red) and type II fibers (blue); (E) Changes in PCr concentration (% initial): the orange line is the simulation result for PCr in the muscle tissue, while black dots with the corresponding line are the experimental data from [92] (mean ± SD for some dots).  Numerical analysis of the model demonstrated the necessity of considering additional transcription factors showing activity 1 to 2 h after exercise for the simulation of genes with a delayed response to exercise (for example, PPARGC1A; [60]). Introducing in the model transcription factor X that is up-regulated immediately after exercise in a CREB-and CRTC-dependent manner allowed us to reproduce the expression of the PPARGC1A gene ( Figure 12K). Our bioinformatics analysis [60] of the transcriptomics data [12], in turn, allowed us to suggest that proteins from the AP-1 family (e.g., FOS and JUN) forming heterodimer complexes with CREB-like transcription factors served as these intermediate regulators (factor X; see details in Supplementary Figure S3 and Table S1). Moreover, an analysis of the transcriptomic [12,80] and ChIP-seq data from the GTRD database [97] revealed that the expression of PPARGC1A via the alternative promoter may be regulated by EGR1 and MYC. Both EGR1 and MYC markedly induced expression 30 to 60 min after an aerobic exercise and had binding motifs in the alternative promoter. Our prediction is supported by experimental data showing that EGR1 expression leads to an increase in PPARGC1A expression in human aortic smooth muscle cells [98,99], while the EGR1 expression promptly and dramatically increased after the stretching of skeletal muscle cells, leading to an increase in the concentration of the EGR1 protein in 3-4 h [100]. On the one hand, MYC positively regulates the expression of all active genes through transcriptional amplification [101][102][103] and chromatin modifications [104,105]. However, an enhancement of its expression negatively impacts PPARGC1A expression [106,107], in particular, in cardiomyocytes [108] and other types of cells where MYC acts as a repressor [109].

The Integrated Modular Model Comprises Three Hierarchical Levels (Metabolic, Signaling, and Gene Expression)
We previously developed a multi-compartmental mathematical model describing the dynamics of intracellular species concentrations and fluxes in human muscle at rest and intracellular metabolic rearrangements in exercising skeletal muscles during aerobic exercise on a cycle ergometer [16]. As an initial model for this study, we used a complex model of energy metabolism in the human skeletal muscle developed by Li and coauthors and considered two types of muscle fibers [32]. We proposed a modular representation of the complex model using the BioUML platform [17]. The modular representation provides the possibility of rapid expansion and modification of the model compartments to account for the complex organization of muscle cells and the limitations of the rate of diffusion of metabolites between intracellular compartments. This feature allowed us to integrate modules of signaling pathways modulating downstream regulatory processes of early response genes and genes with delayed response during exercise and recovery. The validation of the modular model based on a higher number of published experimental data [43,89,90,93,94] (see Supplementary Figure S1) than were used in the original metabolic model [32] showed the validity of the modular modeling approach implemented in BioUML. Furthermore, the integrated modular model provides an absolutely novel in silico platform to predict molecular responses of human skeletal muscle cells to diverse modes of exercise on three hierarchical levels (metabolic, signaling, and gene expression), experimental precise measurements of which are currently methodologically limited or even remain elusive.
In the current state, the model is suitable for testing the plausibility of some physiological hypotheses. For example, the existence of intermediate X factor regulating the expression of the PPARGC1A gene as the example of a delayed response gene in human skeletal muscle has been numerically investigated using different versions of the model: considering direct regulation via the CREB-like factor or taking into account the X factor regulatory role as an intermediate activator of PPARGC1A expression.

Model Constraints and Further Ways for Development
Despite the complexity of the developed modular model, the current version does not consider the influence of many system factors such as hormonal regulation [56,110], the influence of processes in the central nervous system [111,112], feeding mode [113,114], and exercise-induced temperature drift in skeletal muscle [115,116], which hampers the precise quantitative reproduction of abrupt changes at different physiological levels during initial stages of physical exercise. It can be overcome by means of significant modifications on the muscle fiber recruitment model in order to simulate the transient process due to exercise. Some other constraints are described in detail below.
GS activity is regulated through multiple mechanisms, including feedbacks mediated by glycogen, blood glucose concentration, rate of glucose uptake, insulin, epinephrine, and the GS phosphorylation state [46,48,117,118]. However, in the current model, GS activity depends on the glycogen content only. In our model, post-exercise glycogen synthesis is lower than that estimated in the majority of studies [49,119,120] because many factors are omitted, such as feeding and associated rises in blood glucose concentration, rate of glucose uptake, sensitivity to and changes in insulin, etc. At the same time, in our model, glycogen synthesis is higher than that observed during exercise recovery in a fasted state [121]. Additionally, our model does not take into account the effect of muscle fatigue related to the decline in power generated by type II muscle fibres and recruitment of new type II fibres as well as the depletion of muscle glycogen and other substrates. This may play an important role in the simulation of moderate and high intensive and/or long-lasting exercise. Moreover, the focus of this study is related to the recruitment of vastus lateralis muscle fibers and their activation at metabolic, signaling, and gene expression regulation levels as a response to the exercise performed according to a cycle-ergometer or kneeextensor exercises only. These limitations provide a direction for model improvements and should be considered in further works.
Furthermore, the modular nature of the presented model allows the introduction of multiple positive and negative feedbacks between different considered levels: for instance, the impact of kinases altering the activity of enzymes that catalyze reactions of the glycolysis, TCA cycle, and fatty acid oxidation in skeletal muscle [122,123], Ca 2+ -dependent enhancement of glycolytic enzyme activity and mitochondrial respiration [33], and PGC1αdependent regulation of the expression of genes encoding glycolysis and malate-aspartate shuttle enzymes [124].
Our model provides a proof of concept of how dynamic changes at the metabolic level can be linked to gene expression regulation via signaling transduction pathways in skeletal muscles during physical exercises. The modular approach used in the study has demonstrated a methodological basis for qualitative and quantitative development of the complex model including different hierarchical levels of the system organization. The analysis completed during this study allows us to refine the roadmap for further model improvements, linking this in silico version to in vivo skeletal muscle. The roadmap includes an improvement of the motor unit recruitment model, considering the impact of the muscle fatigue on power decrease, and extension of the model by new modules representing system factors, e.g., hormonal regulation and the central nervous system taking into account multiple relationships and feedbacks between different modules of the integrated model.

Conclusions
We developed, for the first time, an integrated model of human skeletal muscle incorporating metabolic, signaling, and gene expression modules. The model enables us to simulate the most important exercise-related signaling (Ca 2+ and AMPK-related signaling) and RNA expression of early response genes (as a result of the activation of transcription factors existing in the cell), as well as the expression of delayed response genes (as a result of the expression of intermediate transcription factors induced immediately after an exercise). The molecular response of skeletal muscle to contractile activity is related to the high number of signaling molecules and genes. The modular nature of the model enables us to add new variables and modules, thereby increasing both the complexity and quality of the model.