Numerical Simulation of a Physiological Mathematical Model of Energy Consumption in a Sarcomere Simulación Numérica de un Modelo Matemático Fisiológico de Consumo de Energía en un Sarcómero

The paradigm of biological systems provides a framework to quantify the behavior of biological processes. Mathematical modeling is one of the analytical tools of biological systems used to reproduce the variables of a system for prediction. This article presents the analysis of muscular contraction, the physiological process responsible of generating force in skeletal muscle, from the point of view of mathematical modeling. The aim is to provide numerical evidences about the force generated by the sarcomere, and the energy required to produce such a force. The proposed scheme includes a model to activate the contractile cycle, based on the action potential that reaches the neuromuscular junction, the calcium release into the sarcoplasm, the contraction response, and the quantification of the energy that the sarcomere requires to perform mechanical work. The results shows that the proposed scheme is acceptable because it reproduces experimental data of force, velocity, and energy reported in the literature. The results of the proposed scheme are encouraging to scale the model at the muscle or muscle group level, in such a way that the quantification of energy can be an alternative to the indirect estimation methods of energy consumption that currently exist.


INTRODUCTION
Biological systems deal with the understanding of biological processes at the systems level. The initial ideas were established by Dr. Norbert Wiener in his book published in the 1940s [1] . However, the execution and confirmation of these ideas did not flourish due to technological (low availability of sensors, actuators, and computer systems), and scientific limitations (for example, the theory of nonlinear systems was in its infancy). It was not until early in this century that these ideas were taken up by the scientific community encouraged, primarily, by advances in molecular biology. The current availability of high-performance computer systems capable of processing copious amounts of data and of information-processing methods like machine learning, as well as the development of sensors capable of measuring biological variables in real time, have all fostered advances in biological systems. In addition, modern systems theory provides a broad platform of methodologies for the analysis, mathematical modeling, and control synthesis of highly-complex processes, including linear, nonlinear, continuous, discontinuous, and interconnected behaviors, to name a few [2] .
The study of biological systems proposes a four-part paradigm for understanding a biological process at the systems level [3] . ( (4) Design method: strategies for the physical implementation of the resulting control scheme according to the first three points. The advantage of this paradigm over the qualitative trial-and-error method used so widely in the biological sciences, is that it provides quantitative descriptions of the process that make it possible to predict the behavior of the properties of interest. Thus, it generates quantitative information, or design parameters, for experimental protocols that help optimize characterization methods and/or design and, more generally, our understanding of the biological process involved. Examples of applications of this paradigm can be consulted in the pioneering papers reported by Kitano [2] [3] . The area of biological systems devoted specifically to human health is called systems medicine, a field that studies physiological processes, pathological conditions, and recommended treatments with the goal of providing quantitative elements to optimize medical treatments [4] .
This article discusses a specific case of analysis of a biological process from the perspective of biological systems: energy consumption in a sarcomere, the basic functional unit of the contraction of skeletal muscle.
The expenditure of energy is defined as the number of calories that people utilize to perform their basic vital functions, and to participate in physical activities [5] .
Total energy expenditure depends on the total energy acquired from food. It is expended, approximately, in the percentages depicted in Figure 1 [6] . in the human body [6] .
Today, energy expenditure is measured by indirect calorimetry, a method that estimates the number of calories the human body consumes while performing a physical activity like walking or running by measuring oxygen consumption (VO 2 ) and carbon dioxide production (VCO 2 ). This allows researchers to calculate the total energy expenditure of an entire human body [6] . Regarding knowledge of the structure of the human body, biological systems conceptualizes this as a system of systems, that can be analyzed on different scales of organization: molecular, cellular, tissue, or organ and, finally, as an integrated organism [4] .
Indirect calorimetry gathers data on energy consumption at the level of the organism; that is, the highest level of analysis of the structure of a biological process, but in certain cases it may be important to determine the functioning of this process on another scale.

Muscular contraction
A skeletal muscle is a tissue that specializes in pro-  [7] . Figure   2 illustrates the organization of a skeletal muscle, from the complete muscle down to the sarcomere. structural protein [7] . The regulator proteins, troponin and tropomyosin, form part of the thin filament, together they make up the troponin-tropomyosin complex. The structure of the sarcomere allows it to contract; that is, to shorten itself by overlapping the thick and thin filaments by through the transformation of energy, from chemical into mechanical [5] [6] .
The contractile cycle is activated when an action potential reaches the neuromuscular union and depolarizes the sarcolemma of the muscle fiber, releasing calcium (Ca 2+ ) into the sarcoplasm. Ca 2+ prepares the thin filament so that the thick filament can bond to it through the association of actin proteins with the heads of the myosin. Actin occupies a site related to myosin that is protected by tropomyosin during relaxation. When Ca 2+ is available in the sarcoplasm, it associates with troponin such that the tropomyosin exposes the sites of the myosin-related actin sites to allow the heads of myosin to bond to that protein.
These unions are called crossbridges. The troponin-tropomyosin complex is recognized as the regulator proteins of muscular contraction due to its function of preparing the thin filament to associate with the thick filament [7] .
The contractile cycle refers to the sequence of events that takes place during the movement of the filaments.
It consists of four stages (see Figure 4). In the first stage (1), myosin becomes charged with energy high. The shortening of the sarcomeres in the myofibrils causes the muscle fiber -and then the complete muscle-to contract [7] .
Muscular contractions are classified as either isotonic or isometric. In the former, the force of contraction developed by the muscle is constant, and the length of the muscle changes. This type of contraction generates movement of the joints and the force required to move loads or objects. In the latter, the force of contraction is insufficient to move a load or object; it only  [7] .
generates sufficient force to sustain it, not move it. In this type of contraction, the muscle does not change its size [7] . Next, the mathematical modeling of muscular contraction is revised.

Mathematical model of contraction
Most of the mathematical models proposed in the literature to emulate the mechanical behavior of skeletal muscles are based on the one posited by Hill et al. [8] ; that is, addressing the mechanical response of muscle at the tissue level [9] [10] . Mathematical modeling of the mechanical response of the muscle, in contrast, is based on the physiological principles of the origin of muscular contraction in the sarcomere of the skeletal muscle. This physiological approach has been used in mathematical modeling of cardiac muscle; the approaches ranges from initial models proposing the mechanical response of a single sarcomere of cardiac muscle [11] , multi-scale mathematical modeling of the heart mechanics [12] [13] to current in silico models used in preclinical trials to assess drugs for cardiac diseases [14] . While our approach sets out from current physio- during the contractile cycle [11] . They called the relaxation stage the phase of ATP hydrolysis. Here, the sarcomere is relaxed, the crossbridges are in a weak conformation (not force-generating), and Ca 2+ is not bonded to troponin. They defined the parameter R (in μM) as the number of units of troponin available in this stage.
The bonding of Ca 2+ to troponin defines the activation stage, when the crossbridges prepare to generate force.
The variable that measures the number of troponin units associated with Ca 2+ is A(t) (in μM). The time variation of A is defined by the next differential equation: where the left-side of the Equation (1)   mere [11] . In the next subsection, the equations to compute the force generated by the sarcomere are presented. Force is related to the velocity of the shortening of the sarcomere, V(t), as well as the troponin units T and U defined in Equation (2) and (3).

Force-velocity relation
Setting out from the assumption that each crossbridge is a pseudoviscous Newtonian element, the force generated by the sarcomere is defined as (Assumption 5 in Landesberg et al. [15] ): and V u (in μm/s) is the maximum velocity of the shortening of the sarcomere [16] .
The sarcomere at rest constitutes the stationary state Landesberg et al. [17] accommodated the terms such that Equation (11) is expressed in the form of Hill's force-velocity ratio equation [18] , where F h is the force in the stationary state: and F m is the force generated by the muscle during isometric contraction:  where SL(0) is the initial length of the sarcomere. In addition, the length of the overlap required to calculate the force generated by the sarcomere, defined in Equation (5), is obtained from the length of the sarcomere, given by the following ratio [19] : where Lφ is the length of the filaments of actin and myosin in a simple overlap during contraction.
Likewise, the force generated by the sarcomere due to isotonic contraction at a certain moment (also called transitory force) is given by the expression [17] : where F 0 =F(0) is the initial condition (at t=0) of the force generated by the sarcomere.

Energy consumed
As mentioned in the description mathematical model of contraction of a sarcomere defined by Equations (1)-(4), the number of crossbridges that pass from the weak to the strong conformation is represented by the variable A(t); that is, the troponin molecule linked to Ca 2+ (t) each one of which corresponds to a crossbridge.
The crossbridges in the strong (force-generating) conformation are represented by the variable T(t). Each crossbridge in the strong conformation requires a unit of ATP hydrolysis and the release of phosphate, in order to pass from the weak to the strong conformation [20] [21] . Hence, the rate of energy consumption is determined by three elements: (1) the variable

A(t);
(2) the transition rate of the crossbridges from the weak to the strong conformation (f ); and (3) the length of the overlap (L s ). This process is modeled by the next differential Equation (15): where E ATP is the free energy released from the hydrolysis of a simple ATP molecule, given by [15] where ρ= 1/g 1 .

Model of Ca 2+ release
Finally, we propose an adaptation of the contraction model defined by Landesberg et al. and represented by Equations (1)-(4) [11] . As it was mentioned before, the For this end, we propose an activation scheme based on a mathematical model to reproduce the input-output response of the dynamics of calcium release in the sarcoplasm. Figure 5 shows this proposal in which the objective is to have a dynamic model based on a transfer function, called model of Ca 2+ release, whose output is the calcium currents through the membrane, this term (I in -I out ). The scheme is based on a transfer function with input data form a train of pulses (AP), the output data correspond to the total Ca 2+ current (I in -I out ) reported by Beuckelmann et al. [22] .
In turn, this term is the one that initiates the contraction cycle defined by the model defined in Equations  The input data of the model of Ca 2+ release is represented by a train of pulses that emulates the AP which reaches the neuromuscular junction. Regarding output data, existing literature has no separate measurements of the currents of Ca 2+ entering (I in ) or leaving (I out ) the sarcoplasm, but experimental studies have measured the total influx of Ca 2+ current, that is (I in -I out ). Thus, the output data used to compute the model of Ca 2+ release are experimental measurements of the total Ca 2+ current (I in -I out ) reported by Beuckelmann et al. [22] .

RESULTS AND DISCUSSION
To devise the model that reproduces the dynamics of total Ca 2+ current in the sarcoplasm, the diagram in  Table 1.        Likewise, the variables A(t), T(t) and U(t) tend to return to their initial values in that period, as Figure 7 shows.     Figure 11 shows the force-length ratio of the shortening of the sarcomere; that is, the muscle's capacity to generate force regardless of the degree of shortening [24] [25] . The amount of tension generated by the muscle depends on how much it can contract or shorten during stimulation.
The orange line represents the result of the simulation by Equations (17) and (18). Figure 11 compares the results to the simulation reported in Tortora et al. [7] , who proposed only minimum and maximum values of the force generated by the sarcomere (image in green). A sarcomere in a relaxed state measures 2-2.2 μm. Figure 11 shows that the maximum force occurs when the sarcomere returns to its normal position (at rest). The normal length range of the sarcomere during the contractile cycle runs from 1.6-2.6 μm. Figure 11 also reveals that the force which decreases the length of the sarcomere is outside the normal range; that is, the length exceeds the value at rest (eccentric contraction), as occurs in extreme contractions. Figure 12 presents the force-velocity ratio. Clearly, during the velocity of shortening (V(t)>0) -that is, contraction-the force tends to remain at its initial or minimum value, before increasing during relaxation.  Finally, the energy consumed by the sarcomere can be seen in Figure 13. The negative results are due to the sign of the velocity of shortening. Here, after 0.5 s, the rate of consumption is small and does not vary, since in this period the Ca 2+ remains in its basal state, and there is no contraction.

CONCLUSIONS
The study of the mechanical properties of the sarcomere of skeletal muscle based on the biological systems paradigm is interesting for characterizing and quantifying the dynamic behavior of muscle cells.
Moreover, the ratio between the mechanical response and the energy consumed by the sarcomere can provide valuable information on the efficient use of energy in cells. The results obtained from the numerical simulations of force are considered acceptable since they coincide with the maximum and minimum ranges of force generated by a sarcomere, according to experimental values reported in the literature [26] [27] .
The availability of experimental data on the maximum and minimum force that the sarcomere can perform is a disadvantage for the validation of the dynamic response of the proposed model. Therefore, to overcome this disadvantage, experimental studies at sarcomere level must be carried out to account with data to validate the dynamical behavior of force generated by the sarcomere. The findings on the length of the sarcomere during the contractile cycle are also deemed acceptable. Experimental reports indicate that the length of a sarcomere at rest phase is 2-4 μm, but when this exceeds the normal state of relaxation, during an extreme contraction, length may decrease to just 1 μm [7] . input are high, the force generated is greater, but that when Ca 2+ is at its basal value, the force remains at its minimum value. Once the model of energy consumption at the level of sarcomere is fully resolved, future work will focus on modeling the coupling rules between sarcomeres that allows the mechanical response of a complete muscle fiber, and on proposing rules for the recruitment of fibers to model the mechanical response and energy consumption of a complete muscle. The challenge seems achievable, since currently studies have been reported addressing this problem of inter-scale modeling, for example, the one reported by Marcucci et al. [28] , where they propose a scaling from muscle fiber to full muscle just considering the mechanical response, without considering energy consumption.
The significance of solving the mathematical modeling of energy consumption by skeletal muscle would be to have quantitative methods, rather than indirect estimations, of a person's energy consumption. This could be of relevant importance in pathologies related to energy management in the human body, from overweight, obesity, metabolic syndrome and diabetes. In addition to the representation in mathematical models of energy consumption, it would allow the design of patient-oriented energy consumption optimization schemes that could be useful, for example, for high-performance athletes.

AUTHOR CONTRIBUTIONS
All authors participate equally in the writing, reviewing, and editing of the present manuscript. K. G. F. R. oversaw the methodology, software, and validation analyses, carried out the investigation and wrote the original draft. D. E. P. G. carried out software analysis and visualization of results. G. Q. C. developed con-ceptualization of the study as well as formal analysis, obtained resources for the study and oversaw the project.

DECLARATION OF INTERESTS
The authors declare that they have no conflict of interests.