INTRODUCTION
ATP is the major energy source for all living organisms and is indispensable for progressing many biological reactions such as muscle contraction, membrane transport, cell migration, and metabolic reactions. ATP demand changes dramatically with cellular activity. For example, in skeletal muscle, ATP is rapidly consumed by muscle contraction during exercise [
1]. In order to avoid ATP depletion during muscle contraction, at the same time, ATP is quickly produced to compensate ATP consumption and to maintain sufficient ATP level. ATP consumption during muscle contraction is mainly caused by actomyosin-ATPase in the sliding motion between actin filaments and myosin filaments following the action potential occurring in the cytoplasm, and Ca
2+-ATPase in control of cytoplasmic Ca
2+ level [
1]. Consumed cytoplasmic ATP is buffered by phosphocreatine (PCr). As a result, PCr becomes creatine in cytoplasm, and mitochondrial ATP is consumed by creatine kinase (CK) localized mitochondria in order to regenerate PCr. ATP production is due to glycolysis in cytoplasm and oxidative phosphorylation (OXPHOS) in mitochondria. Mitochondrial ATP production is triggered by increase in mitochondrial Ca
2+ level as a consequence of increase in cytoplasmic Ca
2+ level by the action potential, which increases the activity of the TCA cycle dehydrogenases and OXPHOS components [
2,
3]. The rate of ATP production and consumption have thus far been indirectly measured. ATP concentration has been also measured during muscle contraction, but these were snapshot data at cell population level [
4]. In order to clarify the quantitative mechanism of ATP production and consumption during muscle contraction, it is necessary to directly measure mitochondrial ATP level and its time course data during muscle contraction from the same individual myotubes through the experiment.
We have developed genetically-encoded Förster resonance energy transfer (FRET)-based ATP indicators that directly measures mitochondrial ATP level [
5,
6]. The ATP indicator, called ATeam, in which the
subunit of F
0F
1-ATP synthase is sandwiched between cyan fluorescent protein (CFP) and yellow fluorescent protein (YFP). Specific binding of ATP to the
subunit in ATeam induces the conformational change of the
subunit from extended to retracted form, and thereby FRET efficiency in ATeam increases. For the targeting ATeam into mitochondrial matrix, duplex of the mitochondrial targeting signal of cytochrome c oxidase subunit VIII was fused to the N terminus of the ATP indicator. The fused indicator, called mitAT1.03, was properly localized to mitochondrial matrix and allowed measurement of mitochondrial ATP level [
5,
7].
In this study, we established C2C12 cells, derived from mouse myoblasts, stably expressing mitAT1.03 ATP probe. By differentiation induction of C2C12 cells, C2C12 cells were differentiated into myotubes, which is an
in vitro cell line model of skeletal muscle functions [
8–
10]. We performed live cell imaging of mitochondrial ATP level in differentiated C2C12 myotubes at a single-cell level with various frequency of electrical stimulation that triggers intracellular Ca
2+ increase and myotube contraction. Based on the time series data of ATP levels obtained by live cell imaging, we constructed a mathematical model for each myotube, and quantitatively analyzed ATP level by balance of ATP production and consumption. We found that differentiated myotubes were clustered into two populations with distinct kinetics of mitochondrial ATP production and consumption, suggesting the possibility that the distinct kinetics reflect the property of fast and slow muscle fibers [
11,
12].
RESULTS
Mitochondrial ATP levels by EPS separated into two clusters, strong response and weak response clusters
The molecular mechanism of mitochondrial ATP production and consumption during muscle contraction has been elucidated [
13] (Fig. 1A), however quantitative control mechanism that maintains ATP level by balance of rapid ATP consumption and production during muscle contraction remains unclear. In this study, we obtained time series data of mitochondrial ATP using mitAT1.03 in individual differentiated C2C12 myotubes by electrical pulse stimulation (EPS) (Fig. 1B). We constructed a mathematical model of mitochondrial ATP for each individual myotubes and analyzed the quantitative control mechanism of mitochondrial ATP production and consumption.
We established a C2C12 cell line stably expressing mitAT1.03, and induced differentiation into myotubes. At first, we confirmed that mitAT1.03 is really localized in the mitochondria by staining differentiated C2C12 myotubes with MitoTracker (Supplementary Fig. S1). Then, the differentiated C2C12 myotubes were stimulated by EPS with 5 different frequencies and FRET ratios of mitAT1.03 were measured to quantify the mitochondrial ATP level (Fig
. 2). The number of stimuli was set to five so as not to damage the myotubes. We developed an automated cell segmentation method of differentiated C2C12 myotubes and used this method for quantification of the FRET signal in individual myotubes (see Materials and methods) [
10]. We performed hierarchical clustering analysis with the obtained time course, which resulted in two clusters (Fig. 2A). Cluster 1 showed strong transient increase by EPS, whereas cluster 2 showed weak transient increase (Fig. 2A, B).
The regulatory pathway of ATP production is faster and more sensitive than that of consumption
Mitochondrial ATP levels are regulated both by ATP production and consumption. For example, strong ATP increase depend on higher production rate and/or lower consumption rate. To examine the contribution of ATP production and consumption to the different ATP response between clusters, we constructed a mathematical model in which mitochondrial ATP levels are regulated by the regulatory pathways of ATP production, Y, and consumption, Z, and these pathways are independently activated by EPS, X (Fig.3A, see Section “Materials and methods” for detail).
We estimated parameters of mathematical model of mitochondrial ATP in individual myotubes (Table 1, Additional information) and simulated time course of mitochondrial ATP level and flux of ATP production and consumption (Fig. 3B, Supplementary Fig. S2).
Each estimated model parameter was compared between ATP production and consumption (Fig. 4A). Among the parameters, the time constant parameters, the apparent dissociation constant parameters, and the basal parameters were significantly different. For the time constant parameter, t, the time constant of regulatory pathway of ATP production, t1, was 0.83 at median, and the time constant of regulatory pathway of ATP consumption parameter, t2, that was 4.89, indicating that regulatory pathway of ATP production is approximately 6 times faster than that of consumption. For the apparent dissociation constant, Kd, the apparent dissociation constant for regulatory pathway of ATP production (Kd1) was 0.99 at median (EC50 was 0.98), and that for regulatory pathway of ATP consumption (Kd2) that was 2.04 (EC50 was 1.61), indicating that regulatory pathway of ATP production is two time more sensitive to the frequency of EPS than that of consumption. Thus, regulatory pathway of mitochondrial ATP production is faster and more sensitive than that of ATP consumption.
Faster ATP production rate resulted in the stronger ATP increase in cluster 1
We also compared each estimated parameter between cluster 1 and cluster 2 (Fig. 4B, Table 1). Among the parameters, k1 of cluster 1 was significantly larger than that of cluster 2 in spite of no significant difference in k2. Because k1 corresponds to a gain in ATP production and k2 corresponds to a gain in ATP consumption, these results indicate that, the strong ATP increase in cluster 1 is because of large gain of ATP production. The time constants, t1 and t2, of cluster 1 were smaller than those of cluster 2. Because t1 and t2 correspond to required time for the regulatory pathways of ATP production and consumption, these results indicate that both regulatory pathways of ATP production and consumption in cluster 1 are faster than that in cluster 2. The Hill coefficient, n2, of cluster 1 was larger than that of cluster 2, indicating that the nonlinearity of regulatory pathway of ATP consumption of cluster 1 is larger than that of cluster 2.
Larger mitochondrial mass in cluster 1 suggests that cluster 1 shows the similar property of slow muscle fiber and cluster 2 shows the similar property of fast muscle fiber
We examined whether the two clusters were separated by cellular properties such as mitochondrial mass and myotube width, and found that mitochondrial mass was larger in cluster 1 than cluster 2 (Fig. 4C). Given that the skeletal muscle in vertebrate consists of slow and fast muscle fibers, and mitochondrial mass of slow muscle fibers is larger than that of fast muscle fibers [
11], this result suggests that cluster 1 shows the similar property of slow muscle fiber, and cluster 2 shows the similar property of fast muscle fiber [
12,
15]. The distinction between slow and fast muscle fibers can be given by specific expression of isoforms of myosin heavy chains (MyHCs): fast fiber expresses fast myosine heavy chain (fMyHCs), and slow fiber expresses slow myosine heavy chain (sMyHCs). Differentiated C2C12 myotubes have been reported that the types of MyHCs expression differs for each myotube [
15]. We tried to examine that the two clusters were derived from the two fiber types, fast and slow muscle fibers, by immunostaining myotubes with antibody specifically recognizing sMyHCs (S58), and antibody specifically recognizing slow and fast MyHCs (F59). However, it is impossible to analyze the relationship directly between the clusters and the fiber types because we cannot measure the time series data of mitochondrial ATP levels and the expression levels of sMyHC simultaneously due to technical limitation. On the other hand, the amount of mitochondria can be measured simultaneously both with mitochondrial ATP levels and the expression levels of sMyHC. Therefore, we tried to examine the relationship between the clusters and expression levels of slow and fast MyHCs via the amount of mitochondria. We stained sMyHCs using antibody specifically recognizing sMyHCs (S58), and found no apparent correlation between the expression level of sMyHCs and YFP intensity (Supplementary Fig. S3). Given that YFP intensity correlates with the mitochondrial mass, the expression level of sMyHCs does not correlate with the mitochondrial mass. We tried to stain both slow and fast MyHCs using antibody specifically recognizing slow and fast MyHCs (F59), but could not obtain the sufficient signal intensity. As a result, we could not obtain the direct evidence that the two clusters correspond to the two fiber types. We examined the correlation between the estimated model parameters and cellular properties, but there was no correlation (Supplementary Fig. S4, 5).
DISCUSSION
In this study, we obtained time series data of mitochondrial ATP level obtained by live cell imaging and constructed mathematical model of mitochondrial ATP level. We characterized model parameters and found that regulatory pathway of ATP production is faster than that of ATP consumption, and that regulatory pathway of ATP production is more sensitive than that of ATP consumption. We performed hierarchical clustering analysis, and found that the time courses of individual myotubes were separated into two clusters showing strong increase (cluster 1) and weak increase (cluster 2). We found the different kinetic characteristics of both ATP production and consumption between cluster 1 and cluster 2.
The detailed mathematical models of ATP production and consumption in mitochondria have been reported [
16–
20]. These models based on biology included many kinds of molecules. These models were at cell population level because they depend on bulk measurement data. On the other hand, our model was simplified, but at single cell level because it depends on time series data obtained by live cell imaging. In models in previous studies, ATP production and consumption are independently activated by stimulation, consistent with our model structure. Although the model parameters in our model cannot be directly compared with those in the previous models, our time series data of mitochondrial ATP obtained by live cell imaging may be useful for model parameter estimation in the previous model.
We found that cluster 1 shows the similar property of slow muscle fiber, and cluster 2 shows the similar property of fast muscle fiber; however, it remains unknown whether both fiber types-like C2C12 myotubes can co-exist in the same dish from a cloned single cell line. Recently, it has been reported that, during the differentiation steps of C2C12 myotubes, sMyHC express in the earlier differentiation step, and fMyHC express in the later differentiation step [
21]. If differentiation of C2C12 myotubes progress with inhomogeneous rate, both fiber types-like C2C12 myotubes can co-exist in the same dish. Therefore, it is possible that cluster 1 corresponds to slow muscle fiber, and cluster 2 corresponds to fast muscle fiber. Slow muscle fibers have the larger oxidative capacity and the smaller glycolytic capacity than fast muscle fibers [
11]. The different kinetics of mitochondrial ATP level between clusters might reflect the different dependence of ATP production on OXPHOS in mitochondria and glycolysis in cytosol between the fiber types.
Interaction between myotubes might affect our conclusion that the kinetics of ATP production and consumption in myotubes by EPS was separated into two clusters. In this study, we ignored the effect of interaction. If interaction affects our conclusion, myotubes in each cluster will show inhomogeneous spatial distribution. Therefore, the spatial cross correlation analysis will suggest whether interaction is the cause of the clusters. However, it is difficult to analyze the effect of interaction because elongated shape of myotube makes it difficult to measure the distances between myotubes for the spatial cross correlation analysis.
We can measure mitochondrial ATP level in skeletal muscle during exercise and construct mathematical model in vivo. It is interesting to compare time course of mitochondrial ATP and model parameters between in vitro and in vivo.
MATERIALS AND METHODS
Establishment of stable cell line
To establish C2C12 cells (kindly provided by Takeaki Ozawa, University of Tokyo, Tokyo, Japan) stably expressing FRET biosensor mitAT1.03 (Plasmid was kindly provided by Hiromi Imamura, Kyoto University, Kyoto, Japan), we used a PiggyBac transposon vector system (System Biosciences, USA). One hundred and fifty µL of Opti-MEM containing 4 µL of Lipofectamine 2000 (Invitrogen, USA) and Opti-MEM containing 1.0 µg of PiggyBac transposon vector clone and 0.2 µg of PiggyBac transposon vector were mixed and incubated for 5 min. Thereafter, 70% confluent C2C12 cells were plated on a 35 mm dish and transfected with the mixture and incubated for 6 h. For selection of infected cells, the cells were cultured with Dulbecco’s Modified Eagle’s Medium (DMEM) (high glucose) with L-glutamine and phenol red (Wako Pure Chemical Industries Limited, Osaka, Japan) containing 10% fetal bovine serum (Nichirei Bioscience Incorporated, Japan) and 20 µg/mL of blasticidin S hydrochloride (Wako, Japan) and cultured until forming the colonies. The colonies were picked and seeded on a Cell Culture Dish 430167. After seeding and proliferation, the cells were stored at a concentration of 1.0×105 cells/mL with bambanker (NIPPON Genetics, Japan). mitAT1.03 is localized in the mitochondria.
Cell culture
C2C12 cells were seeded into 4-well rectangular plates at a density of 3×105 cells/well, with 3 mL of growth medium consisting of DMEM (High glucose) with L-glutamine and phenol red supplemented with 10% fetal bovine serum. This was maintained in an incubator at 37°C under a 5% CO2 atmosphere. On reaching confluence, the medium was changed to a differentiation medium consisting of DMEM (high glucose) with L-glutamine and phenol red supplemented with 2% horse serum (Nichirei Bioscience Incorporated). C2C12 cells were incubated for 7 days in the differentiation medium.
Time-lapse FRET imaging
The medium was changed to Medium 199, Hanks’ Balanced Salts (Life technologies, USA) (M199) with 2% horse serum 14 h before starting observation. Subsequently, differentiated C2C12 myotubes were starved in M199 for 1 h, and then overlaid with mineral oil (Sigma Aldrich) to prevent evaporation of the medium. FRET imaging was performed with a IX83 inverted microscope (Olympus, Tokyo, Japan) equipped with a UPLSAPO10X2 objective lens (Olympus), a ORCAR2 C1060010B CCD camera (Hamamatsu Photonics), a UHGLGPS mercury lamp (Molecular Devices, Sunnyvale, CA) through a 25% neutral density (ND) filter, XF3075 and XF3079 emission filter for CFP and YFP (Omega Optical), respectively, and an MDXY30100TMETA automatically programmable XY stage (Molecular Devices). These images of the myotubes expressing mitAT1.03 were obtained 1-minute intervals. The exposure times were 250 ms for CFP and FRET-YFP images. Myotubes were maintained on a microscope at 37°C. FRETYFP/CFP ratio images were created with MetaMorph software (Universal Imaging, West Chester, PA).
Electrical pulse stimulation
Electrical pulse stimulation was performed using method of Manabe
et al. [
22]. 4-well plate was connected to the electrical stimulation apparatus, a 4-well C-Dish (Ion Optix Corp., Milton, MA, USA), and stimulated by electrical pulses generated by an electrical pulse generator (Uchida Denshi, Hachioji, Japan). Myotubes were stimulated with electrical pulses of 50 V at 0.1, 0.3, 1, 2 or 5 Hz for 10 ms at for given time period during observation with the microscope.
Image analysis
We pre-processed all images to detect individual myotubes using median filter. Individual myotubes were detected by watershed segmentation using a myotube center as the marker. To detect myotubs, time-lapse images of YFP were binarized, and then summed for enhancement of myotubs. All myotubes were detected by binarization of the summed image. To detect the myotube centers, binarized time-lapse images were transformed into distance map images, and then binarized again, finally summed the binarized images. Individual myotube centers were detected by binarization of the summed image. Watershed segmentation was used for detection of individual myotubes using a myotube center as marker.
To quantify intensities of time-lapse images of CFP and YFP, we performed background correction. Because of the region except myotube regions contains undifferentiated cells and the background region, the background intensity in a time-lapse image of CFP and YFP was estimated by two-components Gaussian mixture model. The estimated background intensity was subtracted from the image, then fluorescent intensity of individual myotubes was quantified and calculated FRETYFP/CFP ratio. For more detailed information, see Inoue
et al. [
10].
Subcellular localization and immunocytochemistry
Mitochondria in individual myotubes visualized with YFP of mitAT1.03. or mitochondrial marker, MitoTracker Red FM (Thermo Fisher Scientific, USA).
Myotubes were immunostained on day 7 of differentiation induction [
15]. Myotubes were washed three times with phosphate buffered saline (PBS), pH 7.2, and then fixed with 100% ethanol for 5 min. Myotubes were washed three times with PBS, and blocking solution consisting of 5% horse serum and 2% bovine serum albumin (BSA) in PBS was added to the myotubes for 1 h. S58 and F59 were detected with anti-mouse IgA conjugated to Alexa Fluor 647 and anti-mouse IgG conjugated to Alexa Fluor 568, respectively. The secondary antibodies were diluted 1:100 in PBS containing 10% horse serum, and incubated for 1 h. Myotubes were washed three times with PBS.
Mathematical modeling and numerical simulation
We constructed a mathematical model in which production and consumption of mitochondrial ATP are independently activated by EPS as follows:
where
X represents frequency of EPS.
Y and
Z represent regulatory pathways of mitochondrial ATP production and consumption, respectively.
ATP represents normalized mitochondrial ATP level.
Ybasal and
Zbasal represent the initial values of
Y and
Z, respectively.
Kd1 and
Kd2 represent apparent dissociation constant of Hill equation of
Y and
Z, respectively.
n1 and
n2 represent Hill coefficients of
Y and
Z, respectively. The
nth root of
Kd indicates the half-maximal effective concentration (
EC50).
t1 and
t2 represent time constant of
Y and
Z, respectively.
k1 and
k2 represent rate constant of ATP production and consumption, respectively. The parameters are estimated to reproduce the time course of mitochondrial ATP level during EPS for each myotube (see Section “Parameter estimation”). Numerical simulation was implemented by using MATLAB software (MathWorks).
Parameter estimation
Mitochondrial ATP levels of each myotube were normalized by dividing them by the respective average value of each time course before EPS. The model parameters of each myotube were estimated to reproduce the normalized time course and minimize the objective function value, which is defined as residual sum of square (
RSS) between the normalized time course and the model trajectories.
RSS is given by the equation
where
ATPi and
are the normalized time course and the model trajectories, respectively.
i represents time point of observation with microscopy every minute. The model parameters of each myotube were estimated by 2-step method, composed of global optimization step and local optimization step. In the global optimization step, we used a meta evolutionary programming method, in which we set the number of parents to 500 and the number of generations to 1000, to search the global parameter space for optimized parameters. In local optimization step, we used the nonlinear least squares technique to reach the local minimum from the parameters estimated by the global optimization step [
23]. We regarded this optimized parameter obtained by this 2-step method as the global minimum. We tried this 2-step method 40 times for each myotube and adopted the parameter with the smallest RSS as a parameter of the myotube in the trial. We confirmed that most of 40 parameter sets were similar values and judged that a meta evolutionary programing method has already converged.
Determination of RSS outliers
The outliers of
RSS were detected by the adjusted outlyingness (AO) [
24]. The cutoff value of AO was Q3+ 1.5e3MC ∙ IQR, where Q3, MC, and IQR are the third quartile, medcouple, and interquartile range, respectively. The medcouple is a robust measure of skewness [
25]. The number of directions was set at 8000. Myotubes found to have outlier
RSS (four myotubes) were excluded from further study.
Statistical analysis
The homogeneity of variance of each estimated parameter were confirmed using Levene’s test. The significant differences were confirmed using Brunner-Munzel test. These test results were controlled
FDR using Benjamini-Hochberg procedure [
26]. A
FDR q value of<0.05 was considered statistically significant.
Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature