2.1. Cell Preparation

C2C12 cells (murine myoblasts) were maintained in a growth medium (GM) comprising Dulbecco’s modified Eagle’s medium high glucose (DMEM; Sigma-Aldrich) supplemented with 10% fetal bovine serum (FBS; Biowest) and 1% penicillin-streptomycin (P/S; Sigma-Aldrich). The cells were incubated at 37°C with 5% CO\(_2\). The medium was replaced every two days, and the cells were passaged with trypsin-EDTA (Sigma-Aldrich) upon reaching approximately 80% confluence.

2.2. Formation of Cultured Skeletal Muscle Tissue

Skeletal muscle tissues were fabricated by modifying previously reported methods ^9,22. The culture device comprised a polydimethylsiloxane (PDMS) chamber with longitudinal grooves and anchors at both ends to secure the cultured skeletal muscle tissue. The molds for the chambers and anchors were fabricated using 3D printers (Form 3B\(\texttt{+}\), Formlabs; Mars 5 Ultra 9K, Elegoo; Shape1\(+\)HD, Rayshape) with biocompatible or heat-resistant resins (BioMed Clear, Formlabs; Rigid20, Rayshape). The printed parts were coated with a 2-μm-thick layer of parylene using a chemical vapor deposition system (Parylene Deposition System 2010, Specialty Coating Systems). To fabricate the chamber, a PDMS elastomer and curing agent (Sylgard 184, Dow Corning Toray) were mixed in a \(10:1\) ratio, poured into the molds, and cured at 75°C for 90 min. After release from the mold, the PDMS chamber was obtained. Before use, all the components were sterilized using an ozone/UV sterilizer (CoolCLAVE Blue, Gelantis). To prevent nonspecific tissue adhesion, the groove surfaces were coated with 0.5 wt% 2-methacryloyloxyethyl phosphorylcholine polymer (Lipidure-CM5206, NOF Corporation) in ethanol, heat-treated at 65°C–75°C for 90 min, and subsequently coated with 1 wt% bovine serum albumin (BSA; Sigma-Aldrich). Conversely, the anchors were coated with a 100 μg/mL fibronectin (Sigma-Aldrich) solution in phosphate buffered saline and incubated overnight at 37°C to promote tissue attachment.

For tissue formation, a cell-gel mixture consisting of Matrigel (356231, Corning), GM, fibrinogen (Sigma-Aldrich), thrombin (Sigma-Aldrich), and C2C12 cells (final concentration: 1.5 \(\times\) 10\(^7\) cells/mL) was dispensed into the chamber grooves and around the anchors. After casting, the culture device was incubated at 37°C for 30 min to induce gelation of the cell-gel mixture, followed by the addition of GM supplemented with 1 mg/mL 6-aminocaproic acid (ACA; Sigma-Aldrich). After 1–2 days of culture, the tissues were gently detached from the chamber walls using tweezers. On day 2 of culture, the medium was switched to differentiation medium comprising low-glucose DMEM (Sigma-Aldrich), 2% horse serum (Thermo Fisher Scientific), 1% P/S, and 1 mg/mL ACA. The medium was refreshed every two days throughout the culture period.

2.3. Testing Contractile Motion

On day 9 of culture, contractile motion was evaluated using electrical stimulation. Electrical pulses (1 Hz) were applied to the tissues through gold electrodes using a function generator (WF1947, NF Corporation) and a bipolar amplifier (HSA42011, NF Corporation). During electrical stimulation, tissue motion was recorded using an inverted microscope (IX73, Olympus), and contractile displacement was quantitatively analyzed using motion tracking software (VM-9000, Keyence Corporation). For each sample, the stimulation-synchronized response amplitude was quantified as the peak-to-peak displacement within a fixed response window, and the background noise level was estimated as the standard deviation of the displacement during intervals without stimulation-synchronized responses. In this analysis, a displacement exceeding 20 times the background noise level was used as a quantitative reference for clearly distinguishable contractile motion. The samples labeled as contractile exhibited response amplitudes above this reference level, whereas those labeled as noncontractile showed displacements comparable to the background noise level and no detectable contractile motion.

2.4. Preprocessing of Deep Learning-Based Image Analysis System

Microscopic images of C2C12-based skeletal muscle tissues acquired in a previous study ²² were used to develop a deep learning-based image analysis system. These images were captured once daily using a microscope (IX71N, Olympus) from culture days 1 to 5 (Fig. 1) and utilized for segmentation training, ensemble strategy evaluation, and self-supervised pre-training validation.

To automatically isolate tissue regions, the U\(^2\)-Net model ²³ was employed for tissue segmentation. The training data consisted of 100 images with the ground-truth masks manually annotated using ImageJ, while an additional 20 images were used for validation. All the images were resized to 256 pixels \(\times\) 256 pixels. To compensate for this limited dataset, data augmentation was applied, including horizontal and vertical flipping, random affine transformations (shifts, scaling, and rotations), elastic deformation, grid distortion, brightness and contrast adjustments, and Gaussian blurring.

Binary cross-entropy loss was utilized as the objective function, calculated for each of the multiple-side outputs produced by U\(^2\)-Net, and summed as the final loss. The segmentation performance was evaluated using the Dice similarity coefficient. Optimization was performed via the Adam optimizer with an initial learning rate of \(1\times 10^{-4}\), which was reduced using a ReduceLROnPlateau scheduler (patience \(=5\) epochs). The model was trained for 100 epochs with a batch size of 8. For post-processing, the generated probability maps were binarized, and only the largest connected component was retained to remove small noisy regions. Morphological operations were applied to smooth the boundaries and fill small holes in the segmented regions. After these automated post-processing steps, minor manual refinements were performed only in a small minority of images, roughly on the order of 10%, mainly when anchor regions were partially included in the final segmented masks. These refinements only corrected obvious non-tissue regions, leaving the class labels or downstream classification criteria unaltered.

2.5. Validation of Ensemble Learning Using Culture-Day Identification

Microscopic images from the existing dataset were processed using U\(^2\)-Net for tissue-region extraction and categorized into five classes (day 1 to day 5), reflecting the morphological status, such as tissue compaction and alignment of myotubes. The sample distribution across these classes revealed inherent class imbalances (Table 1), particularly for days 3 and 4, which contained fewer samples than the other stages. To classify these stages, we implemented a soft-voting ensemble comprising three heterogeneous convolutional neural network (CNN) architectures: ConvNeXt-S ²⁴, EfficientNetV2-S ²⁵, and DenseNet201 ²⁶. With the input images resized to 384 pixels \(\times\) 384 pixels, each model was initialized using ImageNet-pre-trained weights and modified with a five-class output layer. These architectures were selected to capture complementary representations based on distinct design principles.

To account for the ordinal nature of culture progression, a hybrid loss function combining soft labels for ordinal regression (SORD)-based KL divergence ²⁷ and Focal Loss ²⁸ was implemented. SORD penalizes distant misclassifications more heavily by generating Gaussian-distributed soft labels (\(\sigma=0.35\)) centered on the ground-truth class. Focal Loss (\(\gamma=2.5\)) was incorporated to prioritize difficult samples by down-weighting the loss from well-classified examples. The final training loss was defined as the weighted sum of these two components, with a 0.7 weight assigned to the SORD-based KL divergence to emphasize the ordinal structure while maintaining the ability to identify the exact ground-truth class through the Focal Loss. Class-specific weights were applied to mitigate the dataset imbalance: day 1 (1.0), day 2 (1.15), day 3 (1.6), day 4 (1.25), and day 5 (1.5). These weights were strategically assigned based on sample size and morphological ambiguity during the intermediate maturation stages to reinforce learning for classes with the highest anticipated difficulty.

For generalization, several data augmentation techniques were employed: random horizontal and vertical flipping with a 50% probability, RandAugment ²⁹ configured to select two operations with a magnitude of 9, and random erasing ³⁰ with 2%–33% of the image area with a 50% probability. For image mixing, MixUp (\(\alpha=0.4\)) ³¹ and CutMix (\(\alpha=1.0\)) ³² were applied with a combined 50% probability, with each method selected equally. Additionally, test-time augmentation (TTA) ³³ was utilized during inference by averaging the classification outputs across original, flipped, and rotated versions of the images.

Optimization was performed using the AdamW optimizer ³⁴ with an initial learning rate of \(7\times 10^{-5}\) and a weight decay of \(1\times 10^{-5}\). A cosine annealing with warmup schedule ³⁵ was applied, where the learning rate was linearly increased from \(1\times 10^{-6}\) up to the initial set value over the first five epochs, followed by a cosine decay back to \(1\times 10^{-6}\) over the remaining epochs. The batch size was set to 16, and early stopping was triggered after 20 epochs, with no improvement in the Macro F1-score. To ensure reproducibility, the random seed was fixed at 42.

A soft-voting ensemble integrated the output probabilities from the three models via a weighted average, with the optimal weights determined by a grid search (step size \(=0.05\)). Evaluation was performed using stratified five-fold cross-validation, and performance was assessed using the Macro F1-score, class-wise F1-scores, and confusion matrix analysis using aggregated out-of-fold classification results (\(n =528\)).

2.6. Self-Supervised Pre-Training Using SimSiam

To compensate for the scarcity of labeled contractile motion data, self-supervised pre-training was conducted using 528 unlabeled culture images preprocessed with U\(^2\)-Net using the SimSiam (Simple Siamese) framework ³⁶. This method learns the domain-specific morphological features by maximizing the similarity between two augmented views of the same image. A stop-gradient operation was employed to prevent the representation collapse without requiring negative samples. EfficientNetV2-S, initialized with ImageNet-pre-trained weights, served as the encoder with input images resized to 256 pixels \(\times\) 256 pixels. The projection and prediction heads were implemented as multilayer perceptrons with an output dimension of 2048, which were specifically configured with three and two layers, respectively. Data augmentation included randomly resized cropping, flipping, 90° rotations, affine transformations (shift, scale), and brightness/contrast adjustments. Optimization was performed using stochastic gradient descent (momentum \(=0.9\), weight decay \(=1\times 10^{-4}\)) with a batch size of 32 for 100 epochs. A cosine decay learning rate schedule was applied, starting at an initial rate of 0.03 after 10 epochs of linear warmup. The resulting pre-trained encoder was utilized as the initial model for the subsequent contractile motion identification task, which was defined as the binary classification of tissues with or without detectable contractile motion under electrical stimulation.

2.7. Identification of Contractile Motion

The presence of contractile motion was evaluated using electrical stimulation and microscopic observations. A dataset of 43 images, comprising 23 positive samples (with contractions) and 20 negative samples (without contractions) was constructed (Fig. 2). The images were preprocessed using the U\(^2\)-Net model to isolate the tissue regions, focusing on the morphological features of the tissue. The task was formulated as a binary classification problem: Class 1 (presence of contractions) and Class 0 (absence of contractions).

To mitigate overfitting risks with small datasets, pre-trained CNNs were employed as fixed-feature extractors rather than being trained end-to-end. Three architectures were evaluated as feature extractors: ConvNeXt-S (input: 384 pixels \(\times\) 384 pixels, 768-dimensional vector), EfficientNetV2-S (input: 384 pixels \(\times\) 384 pixels, output: 1280-dimensional vector), and DenseNet201 (input: 224 pixels \(\times\) 224 pixels, output: 1920-dimensional vector). The feature vectors were obtained via global average pooling ³⁷ after ImageNet normalization. Beyond single models, an ensemble representation was constructed by concatenating the feature vectors from all three architectures, resulting in a 3968-dimensional feature set that captured complementary morphological traits. Additionally, the SimSiam-pre-trained encoder using unlabeled skeletal muscle tissue images was evaluated to determine the efficacy of domain-specific versus ImageNet-based features.

Several classifiers were evaluated, including logistic regression (\(L_1\) and \(L_2\) regularization), support vector machines (linear and RBF kernels), random forest, and gradient boosting. The feature vectors were standardized to zero mean and unit variance before classification. Data diversity was enhanced by applying augmentations including horizontal or vertical flipping, 90° rotations, affine transformations, and brightness adjustments alongside TTA during inference. The dataset was split into 35 training and 8 test samples. Stratified k-fold cross-validation was performed on the training set to identify the optimal classifier and extractor combination, which was retrained on the full training set and evaluated on the held-out test set using the accuracy and F1-score. To ensure statistical robustness given the limited dataset, we evaluated the performance using the average accuracy across 100 independent trials.

2.8. Computational Environment

All models and algorithms were implemented using Python 3.12. Deep-learning models, including U\(^2\)-Net for segmentation and CNN-based feature extractors, were developed using PyTorch. Standard machine learning classifiers such as support vector machines and random forests were implemented using the scikit-learn library. The experiments were conducted on the Google Colab platform using NVIDIA GPU accelerators (e.g., Tesla L4 or A100) to ensure efficient processing and model training.

3.1. Segmentation Performance of U\(^{\boldsymbol{2}}\)-Net

To ensure that the subsequent feature extraction focused exclusively on biological information, a tissue region extraction model was trained using U\(^2\)-Net for 100 epochs. As shown in the loss curves (Fig. 3(a)), both training and validation losses decreased steadily and converged after approximately 50 epochs. The model from epoch 98, which achieved the minimum validation loss of 0.116, was selected as the final model. The optimized model generated high-fidelity probability maps, achieving a Dice coefficient of 0.984 and demonstrating high accuracy in isolating tissue regions.

To evaluate the practical performance of the trained model, the extraction results were visualized using a multi-step refinement process. By applying threshold-based binarization, largest-connected-component extraction, and morphological smoothing to the probability maps generated by U\(^2\)-Net, the tissue regions were precisely isolated (Fig. 3(b)). This visualization confirmed the effective removal of background artifacts, such as the grooves of the culture chamber and anchor connection points, while preserving the tissue structure. The generated masks were applied to the original images to isolate the tissue regions, thereby ensuring that subsequent learning solely utilized morphological features. Furthermore, this automated system improved both efficiency and consistency in data preparation. While manual masking required approximately 5 min per image and displayed subjective variability in defining tissue boundaries, the U\(^2\)-Net-based approach ensured objective and uniform segmentation across all samples. This automated pipeline processed 432 new images, thereby drastically reducing manual labor. In cases where anchor regions were partially misidentified, minor manual refinements were performed to ensure mask quality.

3.2. Ensemble Learning via Culture-Day Identification

To evaluate the effectiveness of the proposed ensemble framework, a validation task was conducted using 528 masked images categorized into five classes depending on the culture day (days 1 to 5). Before ensemble analysis, the performances of three individual architectures, ConvNeXt-S, EfficientNetV2-S, and DenseNet201, were evaluated. A consistent trend across all models was that the highest F1-scores were achieved for day 1, while the lowest were recorded for day 3 (Fig. 3(c)). On day 1, all models exceeded an F1-score of approximately 0.95, reflecting the distinct morphological characteristics at the early stages of culture. Conversely, the F1-scores for day 3 were significantly lower (0.553–0.696), likely due to morphological overlap with consecutive culture days, as well as the limited number of training samples noted in the dataset distribution. Among the individual models, ConvNeXt achieved the best overall performance with a Macro F1-score of 0.796 (Fig. 3(d)).

To determine whether combining diverse architectures could improve the identification stability, a soft-voting ensemble was implemented. The ensemble achieved a Macro F1-score of 0.817 and accuracy of 82.6%, representing improvements of approximately 0.02 and 1.9%, respectively, over the best single model (Fig. 3(d)). Through a grid search (Table 2), the optimal weights were determined to be 0.35 for ConvNeXt, 0.35 for EfficientNetV2, and 0.30 for DenseNet. Notably, although EfficientNetV2 exhibited the lowest individual performance (Macro F1-score \(=0.697\)), it was assigned a weight equal to that of ConvNeXt in the optimal configuration. Ablation analysis showed that excluding ConvNeXt or DenseNet substantially degraded performance (to 0.798 and 0.802, respectively), whereas excluding EfficientNetV2 resulted in a decrease from 0.817 to 0.815. These results suggest that EfficientNetV2 captured complementary image-based morphological representations and distinct misclassification patterns, which enhanced the overall ensemble robustness. Class-wise analysis showed that the ensemble improved the F1-scores across all classes, with the most notable gain on day 4 (from 0.740 to 0.787) (Fig. 3(c)). The confusion matrix revealed that most errors occurred on consecutive days, such as day 2/day 3 (16 and 18 cases) and day 4/day 5 (15 and 12 cases) (Table 3). In comparison, fewer misclassifications were observed between days 3 and 4, suggesting that the morphological changes during this interval were more distinguishable. This may reflect a relatively clearer transition in tissue compaction or overall structural organization around day 4, whereas the boundaries between day 2/day 3 and day 4/day 5 remained more ambiguous. No misclassifications occurred between days 1 and 5, underscoring the clear morphological divergence between the early and late stages.

These results provided two key insights into the proposed system. First, the overall performance gains and optimal weight distributions suggested that heterogeneous architectures captured complementary morphological representations, thereby justifying the soft-voting ensemble approach. These representations may include morphological cues, such as tissue outline, compaction, alignment, and texture patterns associated with myotube-like organization, which are consistent with previously reported structural characteristics of cultured skeletal muscle tissues ²². Second, the concentration of errors between adjacent days (days 2–5) indicated that morphological transitions in the cultured skeletal muscle occurred gradually. These findings validate the efficacy of the proposed feature-extraction framework for the subsequent and more complex task of contractile motion identification.

3.3. Effect of SimSiam Pre-Training

To compensate for the scarcity of labeled data for contractile motion identification, self-supervised pre-training via SimSiam was conducted using unlabeled culture images. This process primarily aimed to enable the encoder to capture the domain-specific morphological features of cultured skeletal muscle tissues without manual annotation. A total of 528 masked images preprocessed using U\(^2\)-Net, which contained structural information related to tissue morphology and culture progression and were therefore suitable for self-supervised representation learning, were used for pre-training. The model was pre-trained for 100 epochs, and the transition of the training loss, defined by the negative cosine similarity, rapidly decreased to approximately \(-0.78\) after 25 epochs during the initial phase (Fig. 4). The loss continued to decline gradually, eventually converging to a final value of \(-0.903\) after approximately 75 epochs. The high degree of similarity between the feature vectors of the two differently augmented views indicates that SimSiam-based pre-training successfully enabled the model to learn invariant representations specific to the morphology of skeletal muscle tissues. By using the pre-trained encoder as the initialization for contractile motion identification, more efficient learning was expected, even with a limited amount of labeled data. These results indicate that SimSiam pre-training converged on the unlabeled skeletal muscle tissue image dataset and enabled the model to learn augmentation-invariant representations. The practical usefulness of the learned representations for contractile motion identification was further evaluated via a subsequent classification task.

3.4. Identification of Contractile Motion

A binary classification task was conducted using labeled microscopic images to identify the presence or absence of contractile motion. To determine the most effective model configuration, the performance of various combinations of feature extractors and classifiers was evaluated (Table 4). The best performance was achieved by combining EfficientNetV2-S with a random forest classifier, which yielded an accuracy and F1-score of 75.0% and 0.800, respectively. In contrast, ConvNeXt-S and DenseNet201 achieved an accuracy of only 62.5%. To evaluate the scalability of the methods validated in the preliminary tasks, the effects of feature-level ensembling and self-supervised pre-training were also analyzed. Unlike the culture-day identification task, in which ensembling improved the performance, the concatenated-feature ensemble achieved only 62.5% accuracy in contractile identification. This result suggests that the high dimensionality of the concatenated vector (3968 dimensions) exceeded the capacity of the model to generalize from the limited training samples, thereby increasing the risk of overfitting. Therefore, future work should optimize feature dimensionality through dimensionality reduction or feature selection under cross-validation, rather than directly using the full concatenated feature vector. Furthermore, the impact of domain-specific pre-training was assessed by comparing the SimSiam-pre-trained encoder against ImageNet-pre-trained models. The SimSiam-based features yielded a test accuracy of only 50.0%, which was 25% lower than that of the ImageNet-pre-trained EfficientNetV2-S (Table 4). This result indicates that while SimSiam learned representations from the 528 unlabeled images, these features were less discriminative for contractile motion than the general-purpose features from ImageNet. The limited diversity and scale of the pre-training dataset possibly hindered the acquisition of task-relevant morphological representations. Overall, increased model complexity via high-dimensional ensembling or small-scale self-supervised learning did not consistently improve the performance under the present data constraints, whereas ImageNet-pre-trained features were more robust.

To examine the classification robustness of the system, an error analysis was performed on the best-performing model from the initial trial (a combination of EfficientNetV2-S and random forest) (Fig. 5). The model correctly classified six of the eight samples. Although 100% accuracy (4/4) was achieved for tissues with contractile motion, only 50% (2/4) of those without motion were correctly identified. Notably, the output probabilities for the correctly classified samples were clustered near the 0.5 decision threshold, ranging from 0.37 to 0.61. This limited confidence suggests that the model may not yet capture sufficiently distinct morphological markers to robustly separate the two functional states. Although a test accuracy of 75.0% was obtained from this single train–test split, an additional evaluation across 100 independent random splits revealed substantially lower mean accuracies (Table 5). The mean accuracies (\(mean{\pm}SD\), \(n =100\)) for EfficientNetV2-S and ConvNeXt-S were \(55.8{\pm}11.6\)% and \(62.1{\pm}10.5\)%, respectively, indicating considerable variability owing to the extremely limited dataset size. A one-sample \(t\)-test against random chance level (50%) confirmed that both results were statistically significant (EfficientNetV2-S: \(p =2.82\times10^{-6}\); ConvNeXt-S: \(p =7.35\times 10^{-20}\)). Nevertheless, the relatively large standard deviations indicate limited stability across different data splits. In addition to the dataset-related limitation, contractile capability cannot be directly determined from unstained static images alone. Although morphological features such as tissue compaction, cellular alignment, and organized myotube-like structures may be associated with contractile function, they do not confirm the functional contraction of the entire tissue. Thus, the present approach should be regarded as a supplementary, noninvasive screening method rather than a replacement for direct contractility assays. Together, these findings indicate that while functional identification from static morphology is feasible at the proof-of-concept level, larger labeled datasets and direct functional validation are required to improve classification reliability. As a practical target, future studies should collect at least several hundred labeled samples across independent culture batches to improve the classification performance and statistical reliability. In addition, future improvements may be achieved through data expansion, transfer learning, incorporation of temporal information from contraction videos, and a more comprehensive performance evaluation using complementary metrics.