A Modeling Insight into Adipose-Derived Stem Cell Myogenesis

Adipose-derived stem cells (ASCs) are clinically important in regenerative medicine as they are relatively easy to obtain, are characterized by low morbidity, and can differentiate into myogenic progenitor cells. Although studies have elucidated the principal markers, PAX7, Desmin, MyoD, and MHC, the underlying mechanisms are not completely understood. This motivates the application of computational methods to facilitate greater understanding of such processes. In the following, we present a multi-stage kinetic model comprising a system of ordinary differential equations (ODEs). We sought to model ASC differentiation using data from a static culture, where no strain is applied, and a dynamic culture, where 10% strain is applied. The coefficients of the equations have been modulated by those experimental data points. To correctly represent the trajectories, various switches and a feedback factor based on total cell number have been introduced to better represent the biology of ASC differentiation. Furthermore, the model has then been applied to predict ASC fate for strains different from those used in the experimental conditions and for times longer than the duration of the experiment. Analysis of the results reveals unique characteristics of ASC myogenesis under dynamic conditions of the applied strain.


Introduction
Adipose-derived stem cells (ASCs) obtained from lipoaspirate tissue provide an easily accessible, abundant source for autologous cells and thus have a great potential to tissue engineering and cell therapies [1][2][3]. ASCs have demonstrated their multi-lineage ability by differentiating into osteogenic, chondrogenic, vascular, neuronal [4,5,3,6] as well as myogenic phenotypes [7]. It has been shown that during myogenesis, ASCs express the same myogenic markers (PAX7/3, Desmin, MyoD, and MHC) and exhibit similar morphological changes (cell alignment and elongation) as satellite cells show in vivo [8,9]. The satellite cells are stem cells that repair and develop skeletal muscle. Despite such promising properties for myogenic purposes, ASCs demonstrate a relatively low differentiation ability; for example, we have shown that these cells subjected to myogenic medium do not express the late marker, MHC [8]. Therefore, there is a need for an improvement in the ASC myogenic capacity. It has recently been shown that the mechanical and biophysical factors, such as cell shape [10], substrate stiffness [11], and surface topography [12,13] play important roles in stem cell fate. Moreover, applied loading (strain) has a substantial effect on stem cell myogenesis as the effects of such strain were explored in the differentiation of ASC and mesenchymal stem cells (MSCs) into smooth muscle cells [14,15]. We have previously considered mechanical cues relevant to the physiological conditions and shown that the application of cyclic uniaxial strains for one hour a day with an amplitude of 10% and frequency of 0.5 Hz can significantly improve ASC myogenesis in vitro [8]. In particular, a significant percent of cells in that study expressed the late myogenic markers, MyoD and MHC, and fused into multinuclear myotubes.
Despite such an outcome, a better understanding and optimization of ASC myogenesis under different conditions are not clear but important. In this regard, a computational model can interpret the experimental observations by using unifying concepts, predict the data beyond the experiment, suggest additional factors to measure, and be extended to study the effects of more complicated culture conditions. Recently, a number of mathematical (computational) models have been developed to predict or explain stem cell fate under different conditions. One approach considers major factors such as a genetic (signaling) network and interprets stem cell differentiation from the mathematical standpoint of the system's bifurcations, i.e. the appearance of additional steady states when an external parameter reaches a critical value [16,17]. A physiological model of hematopoietic cell differentiation [18] has shown that differentiation is governed by the value of a complex parameter that characterizes the ligand/receptor signaling. Another method, suited for kinetics studies, treats stem cell differentiation as a transition through several stages described in terms of fluxes of cell number in a given stage to the next one [19]. The latter approach was previously applied to hematopoietic cells and, among other results, revealed the importance of a feedback signal to make the proposed model more representative of experimental observations [19]. Another important application of kinetic models that incorporate various forms of feedback (signaling) was the analysis of cancer cells [20]. If a model includes multiple interconnected factors, the continuous approach using differential equations may not be effective, and the model of component interaction can be reduced to simpler logical (Boolean) variables [17]. Such findings encourage further implementation of computational methods to ASC differentiation.
We have recently developed an experimental and modeling study to describe in-vitro ASC myogenesis [8]. The proposed model [8] interprets ASC myogenesis as a transition through five stages where each of them is determined by a combination of four myogenic markers (PAX7, Desmin, MyoD, and MHC) whose expression is measured in the experiment. Although the approach [8] reproduced important features of the experimental data, we present here a critical extension of the model that a) broadens its biological framework by incorporating the interactions with and feedback from the cellular environment, b) obtains a better quality approximation of the cell number in each state under static and dynamic conditions, and c) predicts the system's behavior beyond the current experimental conditions, such as for different strains and longer times. Such an advanced model will provide for a better understanding of ASC myogenesis and associated mechanotransduction as well as aid in the design of further experiments to more deeply understand these processes.
In the present paper, we propose a nonlinear model that has first been trained by our experimental data of ASC-myogenesis [8] and then implemented to predict the kinetics for different dynamic conditions. Compared to the model version in [8], we redefine the main stages of the kinetics to accurately reflect the role of the factor PAX7 whose expression first increases to activate the original stem cells and then decreases to promote the process of ASC differentiation. In addition, we assume direct differentiation (without division) in the later stages of the process. We furthermore replace the switches previously related [8] to the particular moments of time with a set of more general time-independent conditions expressed in terms the system variables. Such conditions are biologically associated to signals that arise from the interactions with the myogenic medium and extracellular matrix (ECM). Finally, we introduce a feedback signal that is based on cell density thresholds in the culture. The main features of the model proposed here are shown in introductory Fig 1. We adjust the parameters of our advanced model to reproduce the cell number data from the experimental time course at each stage of the process [8]. We then extend our model to intermediate values of strain and obtain the results for times longer than the duration of the experiment. To do this, we use strain-linear approximations of the model parameters corresponding to 10% (dynamic case) and 0% (static case) strain and compute the kinetics for different strains from 1% through 15%. We found that below the critical value of strain (about 2%) the computed kinetics is similar to the static case. On the other hand, for larger strains, the kinetics exhibit a dynamic pattern where the times of appearance of the later stages, 3 and 4, as well as cell numbers in these stages, are strain-dependent.

Training the Model by Using the Experimental Data
Our extended model is significantly broader than the original version in [8] due to new conditions and parameters included. Thus, we start by adjusting the model parameters using the experimental data from [8, see there also the details of the experiment]. In particular, we use the experimental numbers of cells expressing each of the four myogenic markers (Figure 4c in [8]) as well as the total number of cells (Figure 2b in [8]) reported on days 1, 3, 7, 14, and 21 under static and dynamic (cyclic strain of amplitude 10% and 0.5 Hz frequency) conditions.
Based on the definition of the cell stages in our model (Fig 1), the numbers, n 1 , n 2 , n 3 and n 4 , of cells in stages 1, 2, 3, and 4, respectively, are related to the experimental values n PAX7 , n Des , n MyoD , and n MHC expressing PAX7, Desmin, MyoD, and MHC, respectively, as follows: These relations are broader than those in [8] because they also include the factor PAX7. We determine the left-hand side values in Eqs 1-4 from the experiment and show the results in Fig  2 as open circles (static case) and filled squares (dynamic case). The number of stem cells expressing none of the four factors is not directly reported in the experiment but rather are included in the total cell number, n tot . The modeling data for n 1 , n 2 , n 3 , n 4 , and n tot are presented in Fig 2 via dashed (static cases) and solid (dynamic cases) lines.
After adjusting the model parameters, the results reproduce the major features of the kinetics observed in the experimental data. Indeed, the computed n 1 number undergoes a sharp increase, reaches a maximum at around day 3 of the experiment, and then drops afterwards (Fig 2A). This behavior is typical of both static (open circles) and dynamic (filled squares) cases. The consistent overlap of circles and squares only prior to day 3 indicates that the ASC cultures take different trajectories when strain is applied. Also, the number of cells in the second stage, n 2 , (Fig 2B) monotonically increases throughout the time of the experiment in the static case (dotted line and open circles). In the dynamic case, however, it has a maximum and drops (solid line and filled squares). The number of cells in stage#3, n 3 , stays close to zero until day 7 in the experiment. Similarly, in the model, its deviation from zero is associated with a decreased expression of the factor, PAX7. Previous studies have provided evidence for the mutual inhibition between MyoD and PAX7 [21,22]. Accordingly, with the release of MyoD, PAX7 levels drop, and the number of cells in stage#3 increases, reaches a maximum, and then decreases due to differentiation in the dynamic case (solid line and filled squares in Fig 2C). In the static case, n 3 increases after the release of MyoD too but monotonically increases through day 21 of the experiment (solid line and open circles in Fig 2C). Finally, the cells in the fully differentiated stage#4 are completely absent in the static case ( Fig 2D). In the dynamic case, however, the terminally differentiated cells start appearing at about day 14 of the experiment, their number, n 4 , continue increasing through the final day 21 of the experiment with a slowing rate closer to the end of the experiment(solid line and filled squares in Fig 2D). Such a change in the rate of n 4 increase in the dynamic case is an indication of stabilization for a longer time. Such times are not reached in Fig 2 but they are presented in the predicted data corresponding to a broader dynamic case for different strains below. In the current version of the model, the time at which n 4 deviates from zero depends on whether the factor MyoD (n 3 ) reaches a threshold. Such a threshold is reached in the dynamic case of 10% strain but is not reached in  vivo [23,24] and in vitro [8,9] observations of the final stages of myogenesis. During these stages, myoblasts (in vivo) and precursor cells (in vitro) align, interact, and fuse into myotubes. These processes are likely consequences of the earlier strain application and strain-dependent signaling through the cell-ECM focal adhesions and cell cytoskeleton similar to other adherent cells [25,26]. It is known that MyoD is not the only indicator of the expression of MHC, and a number of characteristic factors (proteins) are involved on the processes cell-cell contact and fusion during myogenesis [27][28][29][30]. Nevertheless, in the experiment with ACS myogenesis [8], the expression of MyoD accompanies cell alignment and fusion, and, for that reason, the MyoD factor has been chosen to represent the precursor compound for the expression of the latest marker MHC. Fig 2E presents the kinetics of the total number of cells in the static (open circles and dashed line) and dynamic (filled squares and solid line) cases, respectively. In both cases, n tot monotonically increases with a slowing rate near the end of the experimental time period. This slowing-down is more pronounced in the dynamic case and is a result of feedback signal limiting the growth in cell numbers when the cell density (total cell number) has reached the critical level. In the static case, such a threshold has not been reached before day 21, but the feedback will have its effect at later times as it is shown below for extended periods of time. The holistic kinetics pictures, including all stages and total numbers of cells, are summarized in Fig 3A and 3B for the static and dynamic cases, respectively.

Static and Dynamic Patterns of Kinetics
The kinetics of cells in the five stages (i = 0,1,2,3,4) are governed by the proliferation, p i , death, d i , and direct differentiation (in the later stages) rates, D i , as well as self-renewal rates, r i . The values of these parameters resulting from the model adjustment to the experimental data are presented in Table 1 (grayed out entries indicate that such parameters are not used). The parameters remain constant and change only at the logical switches as determined by the external signals. In the table, four such switches (absent in the earlier version of the model in [8]) are specified: the moment at which the n 1 -threshold determined by a signal from the myogenic medium (both in the static and dynamic cases); the moment of strain application that triggers the changes in the initial kinetics in the dynamic case governed by the early cell/ECM interaction; the moment when n 1 drops below a threshold that initiates the release of MyoD (both in the static and dynamic cases), and, finally, the time when n 3 reaches a threshold that governs the release of MHC (which is included in both static and dynamic cases but actually implemented only in the dynamic case. The values of the proliferation, p i , death, d i , rates, self-renewal ratios, r i , and direct differentiation rate, D, in five stages of ASC myogenesis under static (no strain) and dynamic (10% strain) conditions.
There are key combinations of the model parameters that determine the pattern of the stage kinetics. In the stages where differentiation occurs via the asymmetric division, such key parameters are as follows: For n 3-kinetics, as r 3 equals one, such a parameter takes the form: Depending on whether the l i -combination is positive or negative, the corresponding n i increases or decreases, respectively. In the beginning, the application of myogenic medium provides a positive l 0 -combination in both the static and dynamic cases; as a result, the original stem cells proliferate and grow in number. This trend is not interrupted in the static case as no strain is applied. Conversely in the dynamic case, the parameters change with strain application and cause l o to switch signs to a negative value, which results in a decrease of the stem cell number, when compared to the static case. The initial increase in the cell number in stage#1 is provided by a significantly positive value of l 1 with almost complete self-renewal and relatively low death rate. After reaching the n 1 -threshold, both the dynamic and static parameters change to those with low self-renewal and higher death rates, resulting in a significantly negative l 1 and a further decrease in n 1 .
According to our concept of the appearance of later myogenic markers, cells do not move into stage#3 until MyoD is released when PAX7 expression (n 1 ) drops to a minimal level. Thus, the equation for n 3 is not integrated, and, consistently, the equation for the previous stage, n 2 , does not include the term responsible for the differentiation into stage#3 (Table 1). After the n 1 -threshold is reached, the equations for n 3 are included in the computations, and the corresponding parameters provide the n 3 -increase, which is more significant in the dynamic case. Similar to the treatment of the MyoD release, the equation for n 4 is not integrated and the differentiation term is not included in the equation for n 3 until n 3 reaches its respective threshold. Since the n 3 -increase is significantly greater in the dynamic case, the conditions of n 4 -release are satisfied at around day 12, after which the number of cells in terminally differentiated stage starts increasing (see more details on the integration of the kinetic equation below).
The main difference between the kinetics of the two cases is that in the static case, the original stem cells continue increasing in number throughout the whole time of the experiment, and the terminally differentiated cells do not appear at all. In contrast, in the dynamic case, the number of stem cells starts decreasing after the strain application, and terminally differentiated cells become dominant toward the end of the experiment, a trend that continues throughout the experiment. Note that the values of the parameters in Table 1 correspond to two particular values, 10% and 0%, of the applied strain. Showing them back-to-back in the table helps to compare the static and dynamic patterns. Importantly, these values are used as the basic points for linear approximations of the model parameters in our study of the strain-dependence of the kinetics of ACS myogenesis.   small strain in Fig 4A, the original stem cells increase in number throughout the whole time interval despite a slight change in the parameters at the moment of strain application. Cells expressing PAX7 (n 1 ) behave similarly to those in the purely static case considered above. Moreover, cells in the late stage#3 increase in number after their release at about day 7, which is the same as in the static case because the n 1 -decrease is the same. The number of cells in stage#3, however, does not reach the threshold required to release the cells in terminally differentiated stage. Thus, by the end of the considered time interval, the fraction of n 3 stem cells becomes the largest one, as there are no cells that reach terminally differentiated stage#4. This is qualitatively similar to the pure static case. Further analysis, not presented, shows that a similar pattern remains up to a strain of about 2%, above which it starts changing. In all cases of strains presented in Fig 4B-4D, the number of cells in stage#0 decreases after strain application, the number of cells in stage#3 reaches its threshold to provide the release of cells in the terminally differentiated stage, and the cell number in stage#4 starts monotonically increasing after that. Moreover, in all the cases in Fig 4B-4D, the number of stem cells becomes small while the fraction of terminally differentiated cells becomes dominant. This is consistent with the dynamic pattern of kinetics. However, analysis reveals differences in the kinetics in these three cases which indicate the strain-dependence of important features of the dynamic pattern. First, the increase in n 3 after its release is steeper, and the moment when n 3 reaches its threshold occurs earlier for larger strains. Such moments happen on days 10 and 19, respectively, in the cases of strain of 15% and 4%. Second, the number of terminally differentiated cells at the end of the same time period is different and higher for larger strains (3.7x10 5 and 2.5x10 5 cells for strains of 15% and 4%, respectively). Also, the rate of increase of n 4 decreases with time in the case of larger strains (compare the results in Fig 4B and 4D). The number of cells in all stages stabilizes with time if we neglect the death rate in stage 4 (not shown). However, the time required until such stabilization is different and is shorter for larger strains. Thus, the main difference between two patterns is the fraction of the remaining stem cells and that of the cell in the terminally differentiated stage. There is a critical value of the strain when the static pattern transforms into the dynamic one, and the dynamic features accordingly become more pronounced for larger strains. Fig 5 summarizes the difference between static and dynamic kinetics of n 0 and n 4 as well as their strain-dependence. Fig 5A shows the straindependence of the fractions of stem cells (n 0 /n tot ) and terminally differentiated cells (n 4 /n tot ) estimated at day 21 (end of the experiment). In the purely static case, the original stem cells take about 50% of the total population, and terminally differentiated cells are absent. With the increase of the applied strain this proportion changes with the fraction of original stem cells quickly diminishing to zero. The terminally differentiated cells appear at about a 3% strain, and ultimately they become the majority and comprise about 80% of the total population. Finally, Fig 5B shows the strain dependence of the time where n o reaches the maximum over the interval of 21 days (thin line) and the time when n 4 starts deviating from zero. For smaller strains (below 2%), n 0 is maximal at the end point of the interval (day 21). Then, for larger strains, when n o starts decreasing after strain application, this time switched to day 3. For strains below 4%, the terminally differentiated cells do not appear by day 21. When a strain of about 4% is applied, the terminally differentiated cells start appearing at an earlier time point, which drops to day 10 for the strain of 15%.

Model and Computational Method
The schematic of our approach is shown in Fig 1, and in this section, we describe the main components of our model and computational details of its implementation. First, we present the system of ODEs for the kinetics of ASC myogenesis, and then, we show how we incorporate the strain dependence into the model. After this, we discuss the computation of the conditions associated with the interaction of the cell with its environment. Then, we describe the computational implementation of the feedback factor due to the limit on cell density. Finally, we discuss the particular behavior of the stage#1 (cell expressing the factor PAX7) in our model.
A flow chart depicting our computational method is presented in Fig 6.

Kinetic Equations
As previously mentioned, our modeling approach is implemented as a system of ODEs whose parameters change at certain moments of time but are constants in between those moments. The stage-to-stage transitions and parameter switching are not fixed in advance but rather obey conditions associated with cell/ECM or cell-cell signaling and are hence expressed in terms of the system variables. Furthermore, the rate of the growth of cell numbers in different stages is constrained due to a physical limit of cell density in the culture, which is implemented via a feedback function of the total cell number. These features make the problem nonlinear, as opposed to its linear counterpart presented in [8].
The equations describing the kinetics of ASC myogenesis are derived as conditions of the flux balance for each of the five stages at a given moment of time. These equations take the form The first terms in the square brackets on the right-hand sides of Eqs 7-9 correspond to the number of self-renewed cells associated with the proliferation process. The second terms in the same brackets correspond to the cell removal due to their death. The third terms in the square brackets in Eqs 8 and 9 correspond to the influx of cells differentiated in the previous stage via cell asymmetric division. The fourth term in the square brackets in Eq 9 and the first terms in Eq 10, respectively, correspond to the direct differentiation from the pre-final to final stage of ASC myogenesis. Finally, the function after the square brackets in Eqs 7-10 corresponds to the feedback factor due to the limit in the cell density.

Strain Dependence
The coefficients in Eqs 7-10 are strain-dependent, and they are assumed to be linear functions of based on the previously determined values for 10% and 0% strain. The values for the parameters corresponding to the 10% and 0% strain as presented above in Table 1 are estimated by matching the experimental data in the dynamic and static cases, respectively. Thus, the coefficients k j (r 0 , p 0 , d 0 , r 1 , p 1 , d 1 , r 2 , p 2 , d 2 , r 3 , p 3 , d 3 , D 3 , and d 4 ) of the system (7-10) are approximated by strain-linear functions as follows conditions and ε is the applied (dimensionless) strain expressed as a decimal. This approximation automatically determines the strain dependence of the parameter combinations, l i , (Eqs 9 and 10) that are key to the pattern of the kinetics of ASC myogenesis. In addition, the conditions that determine the switches between different sets of ODEs are conceptually the same for all strains; however, the actual moments of such switches (except the first of them due to the strain application on day 3) are strain-dependent and they change with the approximation delineated by Eq 11. While the strain approximation of the model rates is monotonic and linear, the resulting kinetics exhibit qualitative changes in the pattern of the system behavior, such as switches from increasing to decreasing in the number of stem cells (Fig 4), switches in moments (days) of reaching the maximal fraction of the stem cell (Fig 5A), switches in moments (days) of appearance of the terminally differentiated cells (Fig 5B), etc. Finally, the linear strain-dependence used here is for a limited range of strain. It is exactly dynamic at 10% and applied up to 15%. Its nonlinear correction as well as additional experimental testing is beyond the scope of the current paper but will be done in the future studies.

Smooth Switches
Step-wise switches in the ODE system parameters may cause the overall solution to have slope discontinuities. In order to avoid this, a smoothening technique is applied to provide a more gradual change in the system parameters. This technique is implemented in four cases associated with strain-application (dynamic case), reaching the n 1 -treshold (static and dynamic cases), n 3 -release (static and dynamic cases), and n 4 -release (dynamic case). The first three of such cases are treated by introducing a gradual change in the parameter combination l i (Eq 5). The last case is treated by using a gradually changing direct differentiation rate, D. Below, we present our analytical approach in the case where the l 1 -parameter changes from a positive value l + to a negative value lnear the moment of strain application (the two other cases associated with a change in parameter l i are treated similarly) The kinetic equation for n 1 is integrated by values of l that depend on what the value of n 1 is. As shown by Eq 12, the progression of l 1 can be described as a piecewise function with four components. Initially, while n 1 is below n 1 Ã (a chosen threshold), "the static" value, l 1 = l + , is used during which n 1 increases. When n 1 reaches n 1 Ã , l 1 begins to linearly decrease as n 1 continues to increase. This behavior continues until n 1 ÃÃ is reached. At n 1 ÃÃ , dn 1 /dt = 0. Moreover, n 1 ÃÃ is slightly above n 1 Ã (1.05-1.1 x n 1 Ã ). Such a relationship between n 1 Ã and n 1 ÃÃ is provided by the choice of the coefficient l in Eq 12. After n 1 ÃÃ , l 1 linearly decreases with a different slope, and n 1 decreases during this stretch as well. After reaching the point where n 1 = n 1 Ã again, the kinetic equation for n 1 is integrated with the "dynamic" value of l 1 = l -, and n 1 continues to decrease. The sketch of the evolution of l 1 and the corresponding behavior of n 1 is shown in Fig  7A and 7B.
In the case of the direct differentiation during the transition of cells from stage#3 to #4, the analytical expression for the gradually changing differentiation rate, D, takes the following form: 8 > > > > > > > > < > > > > > > > > : Once again, the values of D can be described by a piecewise function also with four compartments. Indeed, Eq 13 describes the gradual change of the differentiation rate from 0 (before the n 4 -release) to a positive value, D + (after the direct differentiation starts). A notation consistent with that of Eq 12 has been adopted in the formulation of Eq 13. The changes in the rate D and the corresponding evolution of n 3 are sketched in Fig 7C and 7D.

The Feedback Factor
The feedback factor, f(n tot ) is a function of the total cell number that does not considerably alter cell kinetics while the total cell number is significantly below the threshold, n ‡ . The usage of such a factor prevents unhindered exponential growth, allowing it to better model biological systems. A feedback factor of the following form was chosen to be implemented f ¼ 1=½1 þ e sðn tot Àn z Þ ð 14Þ As the cell count approaches the threshold, however, the feedback factor causes the rate of increase of cell numbers to decrease across all stages. In addition, as soon as the total cell number is slightly above n ‡ , the magnitude of the factor drops to value that is practically zero; in turn, the changes in cell numbers are stabilized for the corresponding stages. In addition, the threshold, n ‡ is also the value that corresponds to a 50% drop in the rate of cell number increase in all five stages. A computational parameter, s, modulates how quickly the rate of cell number increase drops to zero after the threshold cell count is met. In the model, the values of n ‡ % 3-3.5x10 5 and s % 3-4x10 -5 were used. The biological analogue of such a factor can be the result of signals from a protein transiently expressed when the total cell number approaches n ‡ but whose effects persist throughout the myogenic process. The feedback factor presented above is similar to functions used to describe inhibitory-signal processes [17,19].

Stage 1 Suppression
As referenced before, studies have found a mutual inhibition between the PAX7 and MyoD markers [21,22]. Although this biological factor has been incorporated into the model, it was initially found that the n 1 count increases even with MyoD expression in the static and quasistatic cases. It was determined that this can be attributed to the aforementioned feedback factor that provides long-term stabilization. Thus, in order to accurately represent the inhibition, another threshold for was imposed on n 1 , where dn 1 /dt = 0 for any n 1 above the threshold. This was not observed in the dynamic case, and consequently, this additional threshold did not have any effect in the dynamic case and quasi-dynamic cases.

Future Development of the Model
The proposed model describes a fine time-dependent multistage picture of ACS myogenesis, but it assumes a complete spatial homogeneity throughout the area occupied by the cells. To take into account the gradients of the external factors, 3-D geometry, etc., the kinetic equations used here can be incorporated into a PDE model that also accounts for spatial coordinates and the corresponding terms describing the dependence of myogenic processes on these coordinates as it was successfully done in prior cell studies [31]. Alternatively, the concepts of the current method can used for time-discrete rule-based models to compute the time evolving spatial distribution of cells in different stages of differentiation [32]. Another direction of the model development can be a more detailed characterization of cell/ECM and cell/cell signaling obtained based on additional model-stimulated experiments. Finally, a more advanced model can be used to further optimize the experimental procedures, including the most informative set of the reported proteins, time and frequency of the strain application, etc.