A graph neural network model for inferring interindividual variation from experimental biological data

A graph neural network model for inferring interindividual variation from experimental biological data | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article A graph neural network model for inferring interindividual variation from experimental biological data Fuminori Kawano This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7239284/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 12 Nov, 2025 Read the published version in Scientific Reports → Version 1 posted 11 You are reading this latest preprint version Abstract Interindividual variation in biological responses to physiological stimuli is a widely recognized phenomenon. However, effective computational tools for identifying the individual-specific mechanisms remain limited. We present a graph neural network (GNN) model designed to infer hidden molecular and physiological relationships underlying such variation in experimental biological data. To ensure applicability at a laboratory scale, the model was trained on a domain-specific corpus constructed from approximately 65K published studies containing the keyword “skeletal muscle”. The architecture comprises five layers with a multi-head attention mechanism and a multi-layer perceptron, enabling the model to capture both local topological features and directional dominance between connected nodes. The GNN was trained to learn relationships from experimental models to target features, as well as among target features. Using real experimental input consisting of differential gene expression data from mouse skeletal muscle subjected to acute exercise, the model successfully inferred individualized networks, identifying both common and unique paths across individuals based on input experimental context. These results demonstrate the model’s capacity to extract interpretable, individual-specific biological connectivity patterns. The proposed framework serves as a proof of concept for customizable, context-based GNN inference designed to address biological variation at the individual level. Biological sciences/Computational biology and bioinformatics Physical sciences/Mathematics and computing Biological sciences/Neuroscience Biological sciences/Systems biology interindividual variation artificial intelligence machine learning deep learning Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction Deep understanding of interindividual variation, including sex differences, in responses to physiological stimuli such as exercise, nutrition, aging, and pathological conditions is a critical step toward promoting lifelong health and enhancing the quality of human life. The Dunedin Multidisciplinary Health and Development Study, a longitudinal birth cohort study of individuals born in 1972–1973, has revealed substantial variability in the pace of aging among participants 1 , 2 . It further suggests that early-life psychological and nutritional exposures may induce long-term alterations in epigenetic mechanisms that contribute to interindividual differences in aging trajectories. Despite these findings, the processes by which such interindividual variation emerge and develop across the lifespan remain poorly understood. While regular physical activity is widely recognized for its beneficial effects on human health, it is also well documented that considerable variability exists in the outcomes of exercise interventions. For instance, Bamman et al. 3 demonstrated that individuals exhibit heterogeneous responses to resistance training and can be categorized as extreme-, modest-, or non-responders, with these variations not attributable to sex or age. Similar variability has been reported in response to endurance training. Ross et al. 4 reviewed data demonstrating that maximal oxygen consumption increased by only 7% in the lowest responder, whereas the highest responder exhibited an increase of 118% following 24 weeks of treadmill-based endurance training. Furthermore, Bonafiglia et al. 5 quantitatively identified meaningful interindividual differences in trainability by calculating the standard deviation of individual response (SD IR ). They proposed that true variability exists when the SD IR exceeds the smallest worthwhile change, defined as 0.2 times the standard deviation of a control group. Using this criterion, they confirmed significant interindividual variation in the adaptive increases of skeletal muscle citrate synthase activity and capillary density after 4 weeks of cycling exercise training. Collectively, these studies confirm the existence of substantial biological diversity in response to external stimuli. Nevertheless, the underlying mechanisms that give rise to such interindividual variation remain largely elusive. This gap in knowledge may be attributed to the inherent complexity of the biological systems involved, which arises from their multifactorial and person-specific nature, posing a substantial challenge to identifying the determinants of individual responsiveness. Recent advances in artificial intelligence (AI) have enabled in silico analyses to uncover complex patterns and high-dimensional relationships across biological datasets. In bioinformatics, deep learning models trained on genome-wide multi-omics data have been successfully applied to infer the regulatory logic of transcriptional control and to identify key molecular drivers of disease. For example, Xi et al. 6 developed a neural network framework that integrates transcription factor expression with chromatin accessibility signals at putative cis-regulatory elements using a multi-head self-attention mechanism. This model enables accurate prediction of gene expression levels and classification of cell states from single-cell multi-omics data. When applied to type 2 diabetes samples, it further identified candidate regulatory factors with potential therapeutic relevance. Zhang et al. 7 developed a deep clustering framework, scG-cluster, which leverages dual topology-adaptive graph convolutional networks to improve the identification of heterogeneous cell populations in single-cell RNA sequencing data. By jointly modeling feature graphs and cell–cell relationships, their method achieved enhanced accuracy and scalability in predicting fine-grained cellular subtypes across multiple datasets. While these AI-based models provide significant insights into cellular-level variation and regulation, a deeper understanding of interindividual biological variations, particularly in response to physiological stimuli, requires a broader, macro-scale modeling approach. Such a model must integrate molecular networks and physiological parameters with environmental influences that encompass both biological and physiological factors, enabling individualized inference on a mechanistic basis. To investigate the biological mechanisms underlying interindividual variation, it is essential to move beyond static representations of biological systems and consider the dynamic, context-dependent relationships among molecular and physiological factors. While canonical pathway databases such as KEGG have been widely used to interpret high-throughput data by identifying enriched signaling pathways, these static maps often fail to reconcile inconsistencies observed in empirical studies. For instance, experimentally validated changes in gene expression or protein abundance may not align with expected pairwise interactions suggested by public databases. In several cases, only one component of a canonical interaction pair shows significant regulation, whereas the other remains unchanged, suggesting that the activation of biological pathways is conditional and context-specific. Such discrepancies likely stem from heterogeneity in experimental conditions, including differences in tissue types, environmental stimuli, or subject-specific physiological states. To address these limitations, we propose that a dynamic biological network, in which multiple, potentially conflicting factor relationships inferred from diverse experimental contexts are allowed to coexist and interact, is required for accurately modeling individual-specific biological responses. In this study, we present a novel deep learning architecture based on a graph neural network (GNN), designed to infer individualized mechanisms of "bioreaction-variation" — a term that refers to the interindividual variation in biological responses to physiological stimuli. Our GNN framework integrates both molecular-level and physiological-level nodes and is capable of identifying the most plausible mechanistic pathways that explain observed data under specific experimental contexts. Furthermore, the model is also designed to be implementable at a laboratory scale, with a minimized structure that facilitates targeted inference for specific tissues or physiological functions. This design enables context-aware inference of individualized biological networks, thereby offering a scalable and mechanistically interpretable approach for elucidating complex interindividual variation. Methods Overview of bioreaction-variation network The present study introduces a bioreaction-variation network constructed using a GNN framework designed to capture the relationships between experimental models and corresponding physiological or biological parameters, as well as the interactions among these parameters. The GNN architecture comprises two core components: 1) a Model-to-Target interaction layer, which learns the associations between experimental conditions and observed outcomes, and 2) a Target-to-Target interaction layer, which models the interrelationships among the measured parameters themselves. All training data were curated from published experimental studies, with each parameter represented as a change or difference induced by a specific experimental intervention. This structure enables the model, upon input of new experimental data, including detailed conditions and observed parameter changes, to infer relevant biological pathways that mechanistically link experimental conditions and observed physiological responses, enabling individualized interpretation based on patterns learned from previously published studies. Training data collection A total of 65,096 published studies were obtained from PubMed Central. These articles were identified by first searching PubMed with the keyword “skeletal muscle”, and subsequently filtering for those available as free full text in PubMed Central. The main texts of the selected articles were processed using the GPT-4o-mini model to extract relevant experimental contexts. Each study was summarized into a structured format following a fixed JSON schema (Table 1 ). Each JSON entry represented a discrete experimental unit, capturing the details of the experimental model, the physiological or molecular parameters analyzed, and the direction of change. Table 1 Attributes and descriptions for data collection Attribute Description PMID Article ID provided in PubMed Species The species used in the experiment using the common name Age The age or weeks old at the time of the experiment Sex The biological sex of the subject Biosample_main The type of tissue or cells used in the experiment Biosample_detail Additional details about the biosample, such as tissue or cell name Experiment_type The experimental level, e.g., in vivo, in vitro, ex vivo, case study Model_main The main physiological stimulus used in the experiment Model_detail1 The first level of details, such as specific stimulus type Model_detail2 Further details, such as intensity, dosage, duration, or frequency Model_detail3 Additional specific experimental conditions if necessary Timepoint Chronological point at which the sample was collected for analysis Targets Target The specific name of the measured factor Molecule_type Specify the type of molecule, e.g., protein, mRNA Analysis_main The primary measurement method, e.g., western blotting, PCR Analysis_detail Additional details about the measurement method Relation The type of observed change, e.g., increase, decrease Change The degree of change in percentage or fold change Significance The statistical significance using p-values Control The control group used for comparison in "relation" Descriptions shown in this table are excerpts from the actual prompts used in the GPT-4o-mini model for attribute formatting. Note: The attribute was filled by “none”, if applicable answer did not meet in the article. Input data The input data refer to practical datasets used to perform inference with the trained model. In this study, model performance was initially validated using virtual input data constructed to simulate biologically plausible scenarios (see Supplementary data.zip online). Additionally, inference was conducted using real experimental data (see Supplementary data.zip online), derived from RNA sequencing results of individual mouse skeletal muscle samples. These datasets were originally reported in our previous publication 8 and have been deposited on the official journal website associated with that article. To assess interindividual variability, fold changes in gene expression between exercised and non-exercised control mice were calculated across all possible pairwise combinations of individuals. For example, expression profiles of exercise #1 were compared independently to those of control #1, #2, and #3. Genes exhibiting differential responses to exercise were selected based on the following criterion: at least one exercised mouse showed an average fold change in gene expression greater than twofold when compared to all control mice, while the remaining exercised mice exhibited either changes below the twofold threshold or responses in the opposite direction. A total of 15 genes met this criterion and were used for individualized inference in downstream analyses (see Supplementary Fig. S1 online). Embedding To encode experimental context into vector representations, BioBERT (ver. 1.1) 9 was used to tokenize the descriptive elements of each study. This included information related to the experimental models (e.g., species, age, sex, biosample main, biosample detail, experimental type, model main, model detail, and timepoint) as well as target observations (e.g., target, molecular type, analysis main, analysis detail, relation, change, significance, and control). Each set of contextual features was embedded into a 768-dimensional vector representation. To embed model features, particularly those describing experimental species, species-specific biological variation was incorporated through an angular transformation of the original BioBERT embeddings. Specifically, embeddings corresponding to studies involving human, mouse, and rat were adjusted using Rodrigues' rotation formula. This transformation applied a species-specific angular rotation to the original embedding vector, facilitating species-aligned comparisons during inference. Rotation angles were determined based on phylogenetic distances derived from the NCBI Taxonomy database (see Supplementary Fig. S2 online). Specifically, the positions of the 17 species, which together with human, mouse, and rat comprised the top 20 most frequently represented species in the training dataset, were embedded in a two-dimensional space using multidimensional scaling (MDS), and their centroid was treated as the origin. The angular deviation was defined as the angle between the vector from the origin to each of the three major species (human, mouse, and rat) and a reference axis extending from the origin toward the centroid of the 17 species. The radial distance represented the Euclidean distance from the origin to each species point. These values were then used to normalize the original BioBERT embeddings. For target features, embeddings were similarly rotated based on the directionality of reported changes. In cases where the "relation" attribute was labeled as “increase” or “decrease,” the embedding vector was rotated + 90° or − 90°, respectively, using Rodrigues’ rotation formula. This encoding captured the directional semantics of the observed biological changes. Graph construction was performed using PyTorch Geometric (ver. 2.6.1) 10 , in which nodes were defined based on the “model main” categories (model nodes) and individual “target” entries (target nodes) extracted from each publication. Within each study, directed edges were established from the model node to each associated target node. Additionally, undirected edges were formed between all target nodes that were co-analyzed within the same study, resulting in a fully connected subgraph among targets. Each edge was encoded as a 1,536-dimensional tensor by concatenating the 768-dimensional feature embeddings of the source and target nodes. To support node-level learning, node embeddings were also aggregated into mean vectors representing the average features across identical model or target nodes, thereby capturing generalizable node representations across the dataset. Full implementation details, including preprocessing scripts and embedding procedures, are available via the GitHub repository ( https://github.com/fumikawano-lab/Bioreaction-Variation-Network ). GNN model learning GNN training was conducted using a five-layer architecture (Fig. 1 ). Full implementation details, including preprocessing scripts and model training code, are available via the GitHub repository ( https://github.com/fumikawano-lab/Bioreaction-Variation-Network ). To capture Model-to-Target interactions, the first layer employed a multi-head Graph Attention Convolution (GATConv) mechanism 11 , 12 to compute attention weights (attn_W) for each target node. Model and target features, initially encoded as 768-dimensional BioBERT-based embeddings, were linearly transformed into 2,048-dimensional representations composed of 8 attention heads with 256 hidden dimensions per head. Three types of attention weights were computed: attn_m_t: derived from model features and the difference between model and connected target features (i.e., edge features), attn_m_n: similarly based on model-to-target differences, attn_t: derived from target features and edge features. These weights were used to update target node features via a weighted sum: target_x_new = attn_m_t × edge_feature + attn_m_n × target_feature + attn_t × target_feature Subsequent layers (Layers 2–5) implemented a multi-layer perceptron (MLP) to capture Target-to-Target interactions. Each layer comprised fully connected (FC) transformations, batch normalization (BN), and rectified linear unit (ReLU) activation. The second layer standardized all target node features while preserving dimensionality. The third layer reduced the feature dimensionality from 2,048 to 256 through three sequential FC transformations. The fourth layer was designed to model local topological structure via message passing. Message passing was defined as the difference between a given target node’s feature and the average feature of its input neighbors, as derived from the original training data. The mean squared error (MSE) between this reference message passing and the output of the fourth layer was computed as the message passing loss, a component of the total training loss. This reference message passing was then multiplied with each target node’s output feature to weight nodes according to their connectivity and propagation potential. The fifth layer modeled pairwise dominance among target nodes using a domination weight (dom_W). The domination weight was computed as the Euclidean distance between the input feature from a specific source node and the mean feature vector of all input features converging on the same target node, including the input from the source node itself. This distance was then scaled by a sigmoid function. A smaller distance indicates that the input feature closely resembles the aggregate input context of the target node, suggesting a stronger upstream (dominant) influence. In the fifth layer, the updated target feature was obtained through two sequential FC–BN–ReLU transformations, and trained to approximate the reference feature weighted by dom_W. The MSE between the updated features and the reference was used as domination loss, another component of the total training loss. Model optimization was performed using Adam optimizer 13 with a learning rate of 0.001, and backpropagation was carried out over 50 or 200 total epochs. To mitigate gradient bias, total loss was computed as a weighted sum of the two loss components: message passing loss × 2.0 and domination loss × 0.1. Final trained models are publicly available at Google Cloud Storage ( https://storage.googleapis.com/skeletal_muscle/sm_v1/gnn_model/gnn_model_final.pt for 200-epoch model and gnn_model_final_50epoch.pt for 50-epoch model). Network inference To infer individual-specific hidden pathways, network traversal was performed based on user-provided input data comprising experimental model metadata and observed target parameters. These inputs, formatted consistently with the training data structure, were embedded into 768-dimensional feature vectors using BioBERT, following the same preprocessing pipeline used during model training. To initiate the inference process, cosine similarity was computed between the input model vector and all model embeddings in the training dataset. The top five most similar model nodes were selected, and their directly connected target nodes were identified as primary targets, serving as start nodes for subsequent traversal (Step 1 in Fig. 2 ). The aim was to reconstruct viable paths from these start nodes to goal nodes, which corresponded to the target features present in the user-provided input. In Step 2, all possible directed paths between the identified start and goal nodes were enumerated based on the GNN-learned edge connectivity between target nodes. This yielded a candidate path space representing potential mechanistic routes through the network. In Step 3, path refinement was performed by selecting the most contextually relevant edges at each intermediate node. When multiple edges existed between the same node pair (i.e., same source and target nodes but differing feature contexts), the edge whose source feature was most similar to the preceding target node, as determined by cosine similarity, was selected. This ensured feature continuity across the reconstructed path. In Step 4, path quality was assessed using a message passing loss. Initial features at the primary nodes were weighted by the learned attention values (attn_W). For each subsequent node, the predicted feature was updated by multiplying the incoming feature by a sigmoid-scaled Euclidean distance derived from the preceding node, thereby propagating context-specific signal from source to goal. The discrepancy between predicted and input-derived goal node features was computed as the message passing loss for each path. Steps 3 and 4 were embedded in a genetic algorithm framework using the DEAP library 14 . Each candidate path generated in Step 2 was treated as an individual. Fitness was evaluated as a weighted combination of message passing loss and diversity, the latter reflecting the uniqueness of nodes and path topology. Path selection employed the Elitist Non-dominated Sorting Genetic Algorithm (NSGA-II), with mutation applied to 20% of individuals per generation. Mutation involved random substitution of an edge within a path, followed by structural repair using the same edge selection strategy described in Step 2 and 3. If the mutation resulted in higher loss, the original path configuration was retained. Given that the initial path candidates were exhaustively constructed and optimized during Steps 2 and 3, the evolutionary loop was iterated twice, solely to confirm that no further improvement in fitness metrics could be achieved. All source code and execution details for the network inference procedure are available at the associated GitHub repository ( https://github.com/fumikawano-lab/Bioreaction-Variation-Network ). Analysis of individualized networks To identify individual-specific pathways, inferred GNN outputs for each exercised mouse (n = 3) were compared against all non-exercised controls (n = 3). For each exercised mouse, individual networks were constructed by aggregating all predicted paths obtained from pairwise comparisons with each control mouse. A common network was then defined as the set of paths shared across all three individual networks, while unique networks were obtained by subtracting the common network from each individual's network. To assess overall reconstruction accuracy, message passing loss was calculated for each goal node and averaged across the network. In addition, among non-primary and non-goal nodes, the most frequently occurring source node across all inferred networks was identified. For this node, the contribution of each connected target node was assessed by referencing the message passing loss of the specific path (i.e., evolutionary algorithm-derived individual) in which the corresponding edge was included. This provided a direct measure of each target node’s involvement in the reconstruction of input-derived goal node features within individualized networks. Results Overview of training data The training graph was constructed from a total of 27,155 model nodes and 84,723 target nodes. These nodes were interconnected by 383,225 model-to-target edges and 2,475,502 target-to-target edges. Node frequency distributions are summarized in the Supplementary data.xlsx online. Briefly, the most frequently represented experimental models included exercise, high-fat diet feeding, sarcopenia, aging, muscle injury, electrical stimulation, type 2 diabetes, Duchenne muscular dystrophy, diabetes, and isometric contraction. Frequently analyzed target nodes encompassed IL-6, myogenin, MyoD, TNFα, PGC-1α, creatine kinase, insulin, myosin heavy chain, MuRF1, and atrogin1. During model training, the domination loss decreased more rapidly than the message passing loss (see Supplementary Fig. S3 online), a trend that was accounted for by applying differential weighting to these two components in the total loss function. When equal weights were assigned, the message passing loss plateaued early in the training process, limiting further optimization. Although the domination loss approached a minimum around 50 epochs, the message passing loss continued to decline steadily, reducing the difference between the two losses from 0.269 at epoch 50 to 0.228 at epoch 200. These results indicate that training for 200 epochs did not result in overfitting and that both 50-epoch and 200-epoch models were retained for downstream inference to assess potential differences in generalization. Validation of model with virtual input data To evaluate the inference performance of the trained GNN model, two types of virtual input datasets were constructed that differed only in the biological sex of the subjects, with all other experimental conditions held constant. In the networks inferred from the model trained over 50 epochs, key nodes such as myosin heavy chain I, myonuclei, and VO₂ max were frequently identified in both male and female datasets; however, the surrounding edge structures differed noticeably between sexes (Fig. 3 ). Similar nodes were also retrieved from the model trained over 200 epochs, but no clear improvement was observed in the mean distance (average message passing loss) compared to the 50-epoch model, suggesting that inference performance had already stabilized by that point. Network inference with real experimental data To investigate individualized transcriptional responses to exercise, we reanalyzed skeletal muscle RNA sequencing data from exercised and non-exercised mice (n = 3 each) previously reported in our study 8 . Differentially expressed genes between exercised and non-exercised conditions were used to construct input features, which were then inferred through the trained GNN models (trained for either 50 or 200 epochs). Because the contextual input for the experimental model (e.g., exercise protocol, species, tissue) was identical across the three exercised mice, the same primary nodes were consistently identified. However, the inferred networks showed marked interindividual variation in the structure and composition of intermediate and goal nodes (Fig. 4 ). These differences were reflected in the mean message passing loss, which varied across individuals. Notably, inference using the 200-epoch trained model yielded lower average loss values in all individuals compared to the 50-epoch model, indicating improved network reconstruction performance. Given the superior performance of the 200-epoch model, we proceeded to identify common and unique network paths among individuals using this model. A path was defined as a complete sequence of connected nodes from a primary node to a goal node. A total of 27 paths were commonly found across all three individuals, involving key nodes such as AKT, creatine kinase, FOXO3A, mitochondrial morphology, OPA1, and UQCRC2 as intermediate nodes (Fig. 5 A, B). Unique networks were obtained by subtracting this common network from each individual's network (Fig. 5 C). For example, Exer #1 exhibited 65 unique paths, of which 21 overlapped with Exer #2 and 22 with Exer #3 (Fig. 5 A). Exer #2 and Exer #3 shared 35 paths (Fig. 5 A), indicating that, among the three exercised mice analyzed, these two mice shared similar characteristics in the individual-specific regulatory mechanisms governing transcriptional responses to exercise. To further explore key regulatory factors potentially explaining interindividual variation, we focused on intermediate source nodes, excluding primary and goal nodes, and identified those with the highest frequency across all individuals. UQCRC2 emerged as the most recurrent source node in the inferred networks. This high recurrence indicated that UQCRC2 contributed most consistently to distinguishing the hidden network structures across individuals. Figure 6 illustrates the target nodes connected to UQCRC2 in each individual. UQCRC2 formed unique edges to citrate synthase, RPS6, cytochrome c, ANXA2, UQCRC2, CD36, MYH7, MYHC2A, and SIRT3 exclusively in Exer #1. While Exer #2 and Exer #3 shared several downstream nodes, edges from UQCRC2 to triglycerides and free fatty acids were specific to Exer #2, whereas P70S6K, CKMT2, ACTN3, and MYOZ1 were strongly associated with Exer #3. All edges identified from both common and unique paths are listed in the Supplementary data.xlsx online, along with the message passing loss of the paths in which each edge is included. Discussion This study presents a GNN model capable of inferring latent networks that reflect individual-specific responses to physiological stimuli. The model was trained on a dataset constructed from published studies retrieved using the keyword “skeletal muscle”. As of July 2025, approximately 258,000 articles are indexed under this keyword, of which 104,000 are available as free full-text articles via PubMed Central. Our training dataset included 65,096 of these articles, covering more than one-fourth of the accessible literature on skeletal muscle. The predominant experimental models extracted from these publications—such as exercise, high-fat diet, sarcopenia, aging, muscle injury, electrical stimulation, type 2 diabetes, and Duchenne muscular dystrophy—indicate that the dataset reflects biologically relevant contexts closely tied to skeletal muscle research. Accordingly, the extracted parameters and target nodes largely comprised skeletal muscle-related molecules and physiological outcomes. Notably, the model was designed to integrate both molecular and physiological parameters. Among the physiological parameters, non-muscle-related parameters such as VO 2 max, blood glucose, and plasma hormone levels were also included, allowing the reconstruction of networks that reflect both intracellular signaling and systemic physiological outcomes. Rather than focusing solely on inferring molecular cascades, the model learned connections between parameters that frequently co-occur across diverse experimental contexts. This architecture enables the exploration of novel, individualized routes from experimental models to observed phenotypes. This capability is likely supported by the model’s global network structure, which connects parameters across distinct publications through shared nodes and feature similarities. As demonstrated by the inference with virtual input data (Fig. 3 ), even a single biological difference, such as sex, led to the identification of markedly distinct hidden networks. These results highlight the model’s ability to adaptively explore individualized paths. The model’s inference behavior was influenced by the domination weight learned during training, as edge selection was strongly biased by the learned domination-based relationships between nodes. One key difference between the 50- and 200-epoch models appears to lie in the edge directionality and path structure, which became increasingly refined with extended training. Validation using virtual input data demonstrated that the 50-epoch model was sufficient to reconstruct the network, likely due to the simplicity of the input, which included only two goal nodes (Fig. 3 ). In contrast, inference using real biological input data showed improved performance with the 200-epoch model, presumably because the input was more complex, containing unchanged or noisy parameters, thereby requiring deeper learning to refine edge relevance and directionality (Fig. 4 ). These observations collectively suggest that the present GNN functions as a scientifically interpretable, small-scale inference model tailored to skeletal muscle biology, with sufficient capacity to support individualized network inference across a tissue-specific corpus of published studies. Inference using real biological input data derived from the tibialis anterior muscles of mice subjected to a single bout of running following four weeks of training identified Uqcrc2 as a key factor potentially explaining the interindividual variation in transcriptional responses to acute exercise (Fig. 6 ). UQCRC2 (ubiquinol-cytochrome c reductase core protein 2) is a component of mitochondrial respiratory chain complex III, essential for the formation of mitochondrial supercomplexes involving complexes I, III, and IV 15 , 16 . Missense mutations in UQCRC2 have been linked to impaired mitochondrial respiration and inherited human disorders 17 – 19 . Notably, a previous study reported that high volume of high-intensity interval training significantly increased UQCRC2 protein expression in skeletal muscle, contributing to enhanced ATP production through more efficient formation of respiratory complexes 20 . In the individualized networks inferred by the model, UQCRC2 exhibited variable interactions with distinct downstream nodes across mice, suggesting that its local network context contributed differently to each individual’s gene expression response (Fig. 6 ). Importantly, the model's use of domination-based weighting does not imply that UQCRC2 lies upstream in a canonical signaling pathway. Rather, the model hypothesizes directional relationships in which UQCRC2 functions as a dominant influence relative to its connected nodes within the inferred context. For instance, CKMT2, ACTN3, and MYOZ1 were uniquely associated with UQCRC2 in the Exer #3 (Fig. 6 ), implying that these factors might be directly or indirectly modulated by UQCRC2, even though the precise molecular mechanisms linking them remain unknown. This interpretive framework reflects the core capability of the present GNN model: to propose biologically plausible, yet hypothetical, directional associations based on learned patterns across multiple studies. While the inferred relationships should be interpreted as hypothesis-generating rather than confirmatory, the model successfully produced skeletal muscle-specific network structures that align with experimental outcomes and capture individual-level diversity. Taken together, these findings support the utility of this GNN as a domain-focused inference tool with strong contextual validity for skeletal muscle biology. Limitations This study introduces a GNN model designed to infer bioreaction-variation networks within a skeletal muscle–specific corpus. The primary aim was not to propose a universally optimized model, but rather to demonstrate a customizable model architecture that can be implemented and adapted at a laboratory scale to address specific biological questions. The model presented here serves as a proof-of-concept example of such a design. However, the current implementation has several limitations as shown below The model requires input data that represent differential states between experimental and control conditions in order to infer individualized networks. Consequently, datasets comprising only baseline or resting conditions, without any comparative group, are not applicable for inference with this model. Therefore, the extraction of individual-specific features through inference reflects relative interindividual variation within the cohort, including the control subjects used for comparison. Experimental models or analytical parameters for which relevant edges are scarce in the training dataset may lack sufficient connectivity to support robust inference. In such cases, the model may fail to generate appropriate or biologically meaningful network predictions. Because numerical embeddings generated by BioBERT do not directly reflect the magnitude of biological changes in a quantitatively interpretable form, it was necessary to incorporate information on the direction of change, specifically by labeling each relation as either an increase or a decrease. This approach allowed the bioreaction-variation described in the contexts to be translated into vector space. However, capturing such directional changes at the individual level cannot rely on statistical significance testing, which is generally not applicable to single-subject data. In the present study, RNA sequencing data were used, and a twofold change threshold, commonly employed in transcriptomic analyses, was applied to assign directional labels. For other data types or analytical contexts, effect size metrics such as Cohen’s d may provide a more appropriate basis for estimating individual-level deviations from the comparison group. Declarations Competing interests The author declares no competing interests. Ethics declarations Not applicable. Funding The Author received NO FUNDING for this work. Author Contribution FK conceptualized the study, designed the framework of the graph neural network model, wrote all code, validated the model, analyzed the data, drafted the manuscript, and finalized it. Data Availability The datasets generated and/or analysed during the current study are available in the GitHub repository, https://github.com/fumikawano-lab/Bioreaction-Variation-Network. The location of each PyTorch graph dataset on private Google Cloud Storage is also provided in the repository. References Belsky, D. W. et al. Quantification of the pace of biological aging in humans through a blood test, the DunedinPoAm DNA methylation algorithm. eLife 9, (2020). 10.7554/eLife.54870 Belsky, D. W. et al. DunedinPACE, a DNA methylation biomarker of the pace of aging. eLife 11, (2022). 10.7554/eLife.73420 Bamman, M. M., Petrella, J. K., Kim, J. S., Mayhew, D. L. & Cross, J. M. Cluster analysis tests the importance of myogenic gene expression during myofiber hypertrophy in humans. J. Appl. Physiol. 102 , 2232–2239. 10.1152/japplphysiol.00024.2007 (2007). Ross, R. et al. Precision exercise medicine: understanding exercise response variability. Br. J. Sports Med. 53 , 1141–1153. 10.1136/bjsports-2018-100328 (2019). Bonafiglia, J. T. et al. Examining interindividual differences in select muscle and whole-body adaptations to continuous endurance training. Exp. Physiol. 106 , 2168–2176. 10.1113/EP089421 (2021). Xi, X. et al. A mechanism-informed deep neural network enables prioritization of regulators that drive cell state transitions. Nat. Commun. 16 , 1284. 10.1038/s41467-025-56475-9 (2025). Zhang, Y., Feng, X., Wang, Y. & Shi, K. Deep learning powered single-cell clustering framework with enhanced accuracy and stability. Sci. Rep. 15 , 4107. 10.1038/s41598-025-87672-7 (2025). Ohsawa, I. & Kawano, F. Chronic exercise training activates histone turnover in mouse skeletal muscle fibers. FASEB journal: official publication Federation Am. Soc. Experimental Biology . 35 , e21453. 10.1096/fj.202002027RR (2021). Lee, J. et al. BioBERT: a pre-trained biomedical language representation model for biomedical text mining. Bioinformatics 36 , 1234–1240. 10.1093/bioinformatics/btz682 (2020). Fey, M. & Lenssen, J. E. Fast graph representation learning with PyTorch Geometric. ICLR 10.48550/arXiv.1903.02428 (2019). Veličković, P. et al. Graph Atten. Networks ICLR , doi: 10.48550/arXiv.1710.10903 (2018). Zhou, J. et al. Graph neural networks: A review of methods and applications. AI Open. 1 , 57–81. 10.1016/j.aiopen.2021.01.001 (2020). Kingma, D. P., Ba, J. L. & Adam A method for stochastic optimization. ICLR 10.48550/arXiv.1412.6980 (2015). De Rainville, F. M., Fortin, F. A., Gardner, M. A., Parizeau, M. & Gagné, C. DEAP: A Python framework for evolutionary algorithms. Proceedings of the 14th Annual Conference on Genetic and Evolutionary Computation (GECCO ), 85–92, ), 85–92, (2012). 10.1145/2330784.2330799 (2012). Letts, J. A., Fiedorczuk, K. & Sazanov, L. A. The architecture of respiratory supercomplexes. Nature 537 , 644–648. 10.1038/nature19774 (2016). Wu, M., Gu, J., Guo, R., Huang, Y. & Yang, M. Structure of Mammalian Respiratory Supercomplex I(1)III(2)IV(1). Cell 167 (e1510), 1598–1609. 10.1016/j.cell.2016.11.012 (2016). Burska, D. et al. Homozygous missense mutation in UQCRC2 associated with severe encephalomyopathy, mitochondrial complex III assembly defect and activation of mitochondrial protein quality control. Biochim. et Biophys. acta Mol. basis disease . 1867 , 166147. 10.1016/j.bbadis.2021.166147 (2021). Gaignard, P. et al. UQCRC2 mutation in a patient with mitochondrial complex III deficiency causing recurrent liver failure, lactic acidosis and hypoglycemia. J. Hum. Genet. 62 , 729–731. 10.1038/jhg.2017.22 (2017). Miyake, N. et al. Mitochondrial complex III deficiency caused by a homozygous UQCRC2 mutation presenting with neonatal-onset recurrent metabolic decompensation. Hum. Mutat. 34 , 446–452. 10.1002/humu.22257 (2013). Granata, C. et al. High-intensity training induces non-stoichiometric changes in the mitochondrial proteome of human skeletal muscle without reorganisation of respiratory chain content. Nat. Commun. 12 , 7056. 10.1038/s41467-021-27153-3 (2021). Additional Declarations No competing interests reported. Supplementary Files SupplementaryFiguresS13.docx Supplementarydata.xlsx Supplementarydata.zip Cite Share Download PDF Status: Published Journal Publication published 12 Nov, 2025 Read the published version in Scientific Reports → Version 1 posted Editorial decision: Revision requested 13 Aug, 2025 Reviews received at journal 04 Aug, 2025 Reviewers agreed at journal 04 Aug, 2025 Reviews received at journal 04 Aug, 2025 Reviewers agreed at journal 04 Aug, 2025 Reviewers agreed at journal 04 Aug, 2025 Reviewers invited by journal 04 Aug, 2025 Editor assigned by journal 04 Aug, 2025 Editor invited by journal 31 Jul, 2025 Submission checks completed at journal 30 Jul, 2025 First submitted to journal 30 Jul, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7239284","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":496508803,"identity":"c08c54a6-bc03-4abd-b22c-3516b60bf108","order_by":0,"name":"Fuminori Kawano","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA30lEQVRIiWNgGAWjYBACCQhlgyKYQIyWNBDN2ECKlsMoWvADyfYzhp8LKs7XMYgdfv7wZxuDPH8Dw7MH+LRI8+QYS884c1uCQTrNsJm3jcFwxgGGdAN8WuQYcjdI87aBtCQYNjO2MTBugPoMtxb+t5t/87adA2pJ/9gIdJg9QS3SErnbgLYcAGrJMWwAOiyRoBbJGe+/WfOcSZZsk84pnM1zTiJ5xmECfpE4n5Z8m6fCjp9fOn3Dxx9lNrb97T1pD/BpgQM2qBEMDMw8aUTpQAbsx0jWMgpGwSgYBcMaAACbZEBzYk6zwgAAAABJRU5ErkJggg==","orcid":"","institution":"Matsumoto University","correspondingAuthor":true,"prefix":"","firstName":"Fuminori","middleName":"","lastName":"Kawano","suffix":""}],"badges":[],"createdAt":"2025-07-29 05:38:20","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7239284/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7239284/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41598-025-23320-4","type":"published","date":"2025-11-12T15:57:03+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":88518187,"identity":"7474a094-8070-4c58-91a5-18916a39abd4","added_by":"auto","created_at":"2025-08-07 09:21:25","extension":"jpeg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":397497,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eGNN model training overview.\u003c/strong\u003e The graph neural network (GNN) model was trained to capture both Model-to-Target and Target-to-Target interactions through a structured, layer-wise architecture. Layer 1 employed a multi-head Graph Attention Convolution (GATConv) to compute attention weights (attn_W) for each Target node based on associated Model nodes. These attention weights are learnable parameters that allow the model to assign dynamic importance to each model-target interaction. The mechanism operated across 8 attention heads, each with a 256-dimensional hidden space. Here, X denotes the 768-dimensional feature vector of a Target node, aggregated as the mean of all features assigned to the same target entity. Layers 2 to 5 focused on Target-to-Target interactions through a multi-layer perceptron (MLP). Layer 2 applied a fully connected (FC) transformation with batch normalization (BN) and ReLU activation, preserving the feature dimension at 2,048. Layer 3 sequentially encoded the feature vector via three FC-BN-ReLU blocks, progressively reducing the dimensionality from 2,048 to 256. Layer 4 captured local structural dependencies by computing message passing, defined as the difference between the mean feature vector of neighboring nodes (X\u003csub\u003e\u003cem\u003en\u003c/em\u003e\u003c/sub\u003e) and the current node’s feature (X\u003csub\u003e\u003cem\u003ej\u003c/em\u003e\u003c/sub\u003e). This layer also included an FC-BN-ReLU block. Layer 5 inferred pairwise domination between connected Target nodes. The domination weight (dom_W) was calculated as the sigmoid of the feature distance between the mean feature of a target node (X\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e) and the incoming feature vector from another node (e\u003csub\u003e\u003cem\u003ej\u003c/em\u003e\u003c/sub\u003e). Two FC-BN-ReLU blocks were used in this layer as well. All learnable parameters including FC, BN, and attn_W are shown in red.\u003c/p\u003e","description":"","filename":"image1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/ad855a070b45815644ba2ac0.jpeg"},{"id":88516822,"identity":"fc5591b8-7471-4bd4-9eb4-fcf7c0fe15fe","added_by":"auto","created_at":"2025-08-07 09:05:25","extension":"jpeg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":199486,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eNetwork inference overview.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eStep 1: The top five models in the training dataset most similar to the input model (based on model feature similarity) were identified. The target nodes directly connected to these models were then designated as primary nodes, serving as the starting points for network inference.\u003c/p\u003e\n\u003cp\u003eStep 2: Using the trained GNN model, all candidate paths were explored from the primary nodes (start nodes) toward the input target nodes, which served as the goal nodes.\u003c/p\u003e\n\u003cp\u003eStep 3: For each source–target node pair along the paths, the optimal edge was selected from a pool of candidate edges. The edge selection criterion required that the source node of the current edge closely matched the target node of the previous edge in the sequence (i.e., e\u003csub\u003e\u003cem\u003en-1\u003c/em\u003e\u003c/sub\u003e[tgt] ≈ e\u003csub\u003e\u003cem\u003en\u003c/em\u003e\u003c/sub\u003e[src]), thereby ensuring topological continuity toward the goal node.\u003c/p\u003e\n\u003cp\u003eStep 4: Message passing was applied along each candidate path. For the first target node, the attention weight (attn_W) from the connected model node was applied. Subsequent nodes received propagated signals scaled by the sigmoid-transformed Euclidean norm of the previous node’s feature (i.e., 𝜎(∥X\u003csub\u003e\u003cem\u003en\u003c/em\u003e\u003c/sub\u003e∥)). The predicted feature of each goal node (X\u003csub\u003e\u003cem\u003egoal\u003c/em\u003e\u003c/sub\u003e) was obtained at the end of the path. The best combination of paths was selected for each goal node by minimizing the Euclidean distance between the predicted node feature and the actual input feature across all possible path combinations.\u003c/p\u003e","description":"","filename":"image2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/12df476c67c229f26e42d2bf.jpeg"},{"id":88516456,"identity":"2594447e-7f18-4037-9bb1-04f7c3457c7e","added_by":"auto","created_at":"2025-08-07 08:57:25","extension":"jpeg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":305990,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eModel validation using virtual input data.\u003c/strong\u003e Network graphs illustrate all inferred paths based on virtual input data in which the only differing attribute was the biological sex of the subjects. Model training was performed for either 50 or 200 epochs. The message passing loss was calculated for the best combination of inferred paths associated with each goal node, and these values were averaged to produce the mean distance (mean dist.).\u003c/p\u003e","description":"","filename":"image3.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/dcf674db0a640faf46e4f462.jpeg"},{"id":88517742,"identity":"148f9a40-a4f6-4361-8848-d3d00f863241","added_by":"auto","created_at":"2025-08-07 09:13:25","extension":"jpeg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":817078,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eNetwork inference using real experimental data.\u003c/strong\u003e Network graphs illustrate all inferred paths for each exercised mouse (Exer #1–3), based on reanalyzed RNA sequencing data from individual tibialis anterior muscle samples obtained in our previous study \u003ca href=\"#_ENREF_8\" title=\"Ohsawa, 2021 #651\"\u003e\u003csup\u003e8\u003c/sup\u003e\u003c/a\u003e. A single bout of treadmill running was conducted after a 4-week training protocol. Gene expression profiles were compared between three exercised mice and three non-exercised controls. See the Methods section for details on the individual comparisons. Model training was performed for either 50 or 200 epochs. The message passing loss was calculated for the best combination of inferred paths associated with each goal node, and these values were averaged to produce the mean distance (mean dist.).\u003c/p\u003e","description":"","filename":"image4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/9c5f1b768c46c512d7df8e94.jpeg"},{"id":88516452,"identity":"6b3809d6-4bb7-415f-893c-42e3d94672e6","added_by":"auto","created_at":"2025-08-07 08:57:25","extension":"jpeg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":376933,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eExtraction of common and individual-specific networks.\u003c/strong\u003e (a) Venn diagram showing the number of shared and unique inferred paths among the three exercised mice (Exer #1–3). (b) The common network represents the overlapping paths found in all three individuals. (c) Unique networks for each mouse were extracted from the total inferred networks presented in Fig. 4 by identifying individual-specific paths.\u003c/p\u003e","description":"","filename":"image5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/390fc4db93b280e375539b1e.jpeg"},{"id":88516458,"identity":"90852af0-5ffc-4e76-8a28-3a2b0fa8d020","added_by":"auto","created_at":"2025-08-07 08:57:25","extension":"jpeg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":117194,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eVariation of UQCRC2-associated edges.\u003c/strong\u003e Heat map illustrating the target nodes connected to UQCRC2 across the three exercised mice (Exer #1–3). Targets were re-ordered by hierarchical clustering. Color intensity reflects the message passing loss of the path including the corresponding edge, with darker blocks indicating lower loss.\u003c/p\u003e","description":"","filename":"image6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/1b26e6083d75b2c882b93633.jpeg"},{"id":96105020,"identity":"8974b362-ee2d-4f59-bfef-00e1233f486e","added_by":"auto","created_at":"2025-11-17 16:07:15","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2795992,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/9a7a391c-781c-4992-872e-e8e3b989bc11.pdf"},{"id":88516827,"identity":"c7d873d2-47cc-4fdd-8041-1af237b076fa","added_by":"auto","created_at":"2025-08-07 09:05:25","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":483362,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFiguresS13.docx","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/71472406e31a1fb27f6c646a.docx"},{"id":88518190,"identity":"566c1d9c-be66-462b-935e-4a092ca11cd6","added_by":"auto","created_at":"2025-08-07 09:21:26","extension":"xlsx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":2321084,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarydata.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/8ea67fd12b39013b92568105.xlsx"},{"id":88516825,"identity":"07613c07-75c0-4ee2-89a9-bacf055129ef","added_by":"auto","created_at":"2025-08-07 09:05:25","extension":"zip","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":7133,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarydata.zip","url":"https://assets-eu.researchsquare.com/files/rs-7239284/v1/eb66e35263c4b79ee4ee31b9.zip"}],"financialInterests":"No competing interests reported.","formattedTitle":"A graph neural network model for inferring interindividual variation from experimental biological data","fulltext":[{"header":"Introduction","content":"\u003cp\u003eDeep understanding of interindividual variation, including sex differences, in responses to physiological stimuli such as exercise, nutrition, aging, and pathological conditions is a critical step toward promoting lifelong health and enhancing the quality of human life. The Dunedin Multidisciplinary Health and Development Study, a longitudinal birth cohort study of individuals born in 1972–1973, has revealed substantial variability in the pace of aging among participants \u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. It further suggests that early-life psychological and nutritional exposures may induce long-term alterations in epigenetic mechanisms that contribute to interindividual differences in aging trajectories. Despite these findings, the processes by which such interindividual variation emerge and develop across the lifespan remain poorly understood. While regular physical activity is widely recognized for its beneficial effects on human health, it is also well documented that considerable variability exists in the outcomes of exercise interventions. For instance, Bamman et al. \u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e demonstrated that individuals exhibit heterogeneous responses to resistance training and can be categorized as extreme-, modest-, or non-responders, with these variations not attributable to sex or age. Similar variability has been reported in response to endurance training. Ross et al. \u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u003c/sup\u003e reviewed data demonstrating that maximal oxygen consumption increased by only 7% in the lowest responder, whereas the highest responder exhibited an increase of 118% following 24 weeks of treadmill-based endurance training. Furthermore, Bonafiglia et al. \u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e quantitatively identified meaningful interindividual differences in trainability by calculating the standard deviation of individual response (SD\u003csub\u003eIR\u003c/sub\u003e). They proposed that true variability exists when the SD\u003csub\u003eIR\u003c/sub\u003e exceeds the smallest worthwhile change, defined as 0.2 times the standard deviation of a control group. Using this criterion, they confirmed significant interindividual variation in the adaptive increases of skeletal muscle citrate synthase activity and capillary density after 4 weeks of cycling exercise training. Collectively, these studies confirm the existence of substantial biological diversity in response to external stimuli. Nevertheless, the underlying mechanisms that give rise to such interindividual variation remain largely elusive. This gap in knowledge may be attributed to the inherent complexity of the biological systems involved, which arises from their multifactorial and person-specific nature, posing a substantial challenge to identifying the determinants of individual responsiveness.\u003c/p\u003e\u003cp\u003eRecent advances in artificial intelligence (AI) have enabled in silico analyses to uncover complex patterns and high-dimensional relationships across biological datasets. In bioinformatics, deep learning models trained on genome-wide multi-omics data have been successfully applied to infer the regulatory logic of transcriptional control and to identify key molecular drivers of disease. For example, Xi et al. \u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e developed a neural network framework that integrates transcription factor expression with chromatin accessibility signals at putative cis-regulatory elements using a multi-head self-attention mechanism. This model enables accurate prediction of gene expression levels and classification of cell states from single-cell multi-omics data. When applied to type 2 diabetes samples, it further identified candidate regulatory factors with potential therapeutic relevance. Zhang et al. \u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e developed a deep clustering framework, scG-cluster, which leverages dual topology-adaptive graph convolutional networks to improve the identification of heterogeneous cell populations in single-cell RNA sequencing data. By jointly modeling feature graphs and cell–cell relationships, their method achieved enhanced accuracy and scalability in predicting fine-grained cellular subtypes across multiple datasets. While these AI-based models provide significant insights into cellular-level variation and regulation, a deeper understanding of interindividual biological variations, particularly in response to physiological stimuli, requires a broader, macro-scale modeling approach. Such a model must integrate molecular networks and physiological parameters with environmental influences that encompass both biological and physiological factors, enabling individualized inference on a mechanistic basis.\u003c/p\u003e\u003cp\u003eTo investigate the biological mechanisms underlying interindividual variation, it is essential to move beyond static representations of biological systems and consider the dynamic, context-dependent relationships among molecular and physiological factors. While canonical pathway databases such as KEGG have been widely used to interpret high-throughput data by identifying enriched signaling pathways, these static maps often fail to reconcile inconsistencies observed in empirical studies. For instance, experimentally validated changes in gene expression or protein abundance may not align with expected pairwise interactions suggested by public databases. In several cases, only one component of a canonical interaction pair shows significant regulation, whereas the other remains unchanged, suggesting that the activation of biological pathways is conditional and context-specific. Such discrepancies likely stem from heterogeneity in experimental conditions, including differences in tissue types, environmental stimuli, or subject-specific physiological states. To address these limitations, we propose that a dynamic biological network, in which multiple, potentially conflicting factor relationships inferred from diverse experimental contexts are allowed to coexist and interact, is required for accurately modeling individual-specific biological responses. In this study, we present a novel deep learning architecture based on a graph neural network (GNN), designed to infer individualized mechanisms of \"bioreaction-variation\" — a term that refers to the interindividual variation in biological responses to physiological stimuli. Our GNN framework integrates both molecular-level and physiological-level nodes and is capable of identifying the most plausible mechanistic pathways that explain observed data under specific experimental contexts. Furthermore, the model is also designed to be implementable at a laboratory scale, with a minimized structure that facilitates targeted inference for specific tissues or physiological functions. This design enables context-aware inference of individualized biological networks, thereby offering a scalable and mechanistically interpretable approach for elucidating complex interindividual variation.\u003c/p\u003e"},{"header":"Methods","content":"\u003cp\u003e\u003cem\u003eOverview of bioreaction-variation network\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe present study introduces a bioreaction-variation network constructed using a GNN framework designed to capture the relationships between experimental models and corresponding physiological or biological parameters, as well as the interactions among these parameters. The GNN architecture comprises two core components: 1) a Model-to-Target interaction layer, which learns the associations between experimental conditions and observed outcomes, and 2) a Target-to-Target interaction layer, which models the interrelationships among the measured parameters themselves. All training data were curated from published experimental studies, with each parameter represented as a change or difference induced by a specific experimental intervention. This structure enables the model, upon input of new experimental data, including detailed conditions and observed parameter changes, to infer relevant biological pathways that mechanistically link experimental conditions and observed physiological responses, enabling individualized interpretation based on patterns learned from previously published studies.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eTraining data collection\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eA total of 65,096 published studies were obtained from PubMed Central. These articles were identified by first searching PubMed with the keyword \u0026ldquo;skeletal muscle\u0026rdquo;, and subsequently filtering for those available as free full text in PubMed Central. The main texts of the selected articles were processed using the GPT-4o-mini model to extract relevant experimental contexts. Each study was summarized into a structured format following a fixed JSON schema (Table\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e). Each JSON entry represented a discrete experimental unit, capturing the details of the experimental model, the physiological or molecular parameters analyzed, and the direction of change.\u003c/p\u003e\n\u003cdiv class=\"gridtable\"\u003e\n\u003ctable id=\"Tab1\" border=\"1\"\u003e\u003ccaption\u003e\n\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\n\u003cdiv class=\"CaptionContent\"\u003e\n\u003cp\u003eAttributes and descriptions for data collection\u003c/p\u003e\n\u003c/div\u003e\n\u003c/caption\u003e\u003ccolgroup\u003e\u003c/colgroup\u003e\n\u003cthead\u003e\n\u003ctr\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eAttribute\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eDescription\u003c/p\u003e\n\u003c/th\u003e\n\u003c/tr\u003e\n\u003c/thead\u003e\n\u003ctbody\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003ePMID\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eArticle ID provided in PubMed\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eSpecies\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe species used in the experiment using the common name\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAge\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe age or weeks old at the time of the experiment\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eSex\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe biological sex of the subject\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eBiosample_main\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe type of tissue or cells used in the experiment\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eBiosample_detail\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAdditional details about the biosample, such as tissue or cell name\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eExperiment_type\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe experimental level, e.g., in vivo, in vitro, ex vivo, case study\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eModel_main\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe main physiological stimulus used in the experiment\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eModel_detail1\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe first level of details, such as specific stimulus type\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eModel_detail2\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eFurther details, such as intensity, dosage, duration, or frequency\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eModel_detail3\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAdditional specific experimental conditions if necessary\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eTimepoint\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eChronological point at which the sample was collected for analysis\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eTargets\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eTarget\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe specific name of the measured factor\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eMolecule_type\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eSpecify the type of molecule, e.g., protein, mRNA\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAnalysis_main\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe primary measurement method, e.g., western blotting, PCR\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAnalysis_detail\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eAdditional details about the measurement method\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eRelation\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe type of observed change, e.g., increase, decrease\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eChange\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe degree of change in percentage or fold change\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eSignificance\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe statistical significance using p-values\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eControl\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eThe control group used for comparison in \"relation\"\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003c/tbody\u003e\n\u003ctfoot\u003e\n\u003ctr\u003e\n\u003ctd colspan=\"2\"\u003eDescriptions shown in this table are excerpts from the actual prompts used in the GPT-4o-mini model for attribute formatting. Note: The attribute was filled by \u0026ldquo;none\u0026rdquo;, if applicable answer did not meet in the article.\u003c/td\u003e\n\u003c/tr\u003e\n\u003c/tfoot\u003e\n\u003c/table\u003e\n\u003c/div\u003e\n\u003cp\u003e\u003cem\u003eInput data\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe input data refer to practical datasets used to perform inference with the trained model. In this study, model performance was initially validated using virtual input data constructed to simulate biologically plausible scenarios (see Supplementary data.zip online). Additionally, inference was conducted using real experimental data (see Supplementary data.zip online), derived from RNA sequencing results of individual mouse skeletal muscle samples. These datasets were originally reported in our previous publication \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e and have been deposited on the official journal website associated with that article. To assess interindividual variability, fold changes in gene expression between exercised and non-exercised control mice were calculated across all possible pairwise combinations of individuals. For example, expression profiles of exercise #1 were compared independently to those of control #1, #2, and #3. Genes exhibiting differential responses to exercise were selected based on the following criterion: at least one exercised mouse showed an average fold change in gene expression greater than twofold when compared to all control mice, while the remaining exercised mice exhibited either changes below the twofold threshold or responses in the opposite direction. A total of 15 genes met this criterion and were used for individualized inference in downstream analyses (see Supplementary Fig. \u003cspan class=\"InternalRef\"\u003eS1\u003c/span\u003e online).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eEmbedding\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo encode experimental context into vector representations, BioBERT (ver. 1.1) \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e was used to tokenize the descriptive elements of each study. This included information related to the experimental models (e.g., species, age, sex, biosample main, biosample detail, experimental type, model main, model detail, and timepoint) as well as target observations (e.g., target, molecular type, analysis main, analysis detail, relation, change, significance, and control). Each set of contextual features was embedded into a 768-dimensional vector representation.\u003c/p\u003e\n\u003cp\u003eTo embed model features, particularly those describing experimental species, species-specific biological variation was incorporated through an angular transformation of the original BioBERT embeddings. Specifically, embeddings corresponding to studies involving human, mouse, and rat were adjusted using Rodrigues' rotation formula. This transformation applied a species-specific angular rotation to the original embedding vector, facilitating species-aligned comparisons during inference. Rotation angles were determined based on phylogenetic distances derived from the NCBI Taxonomy database (see Supplementary Fig. \u003cspan class=\"InternalRef\"\u003eS2\u003c/span\u003e online). Specifically, the positions of the 17 species, which together with human, mouse, and rat comprised the top 20 most frequently represented species in the training dataset, were embedded in a two-dimensional space using multidimensional scaling (MDS), and their centroid was treated as the origin. The angular deviation was defined as the angle between the vector from the origin to each of the three major species (human, mouse, and rat) and a reference axis extending from the origin toward the centroid of the 17 species. The radial distance represented the Euclidean distance from the origin to each species point. These values were then used to normalize the original BioBERT embeddings.\u003c/p\u003e\n\u003cp\u003eFor target features, embeddings were similarly rotated based on the directionality of reported changes. In cases where the \"relation\" attribute was labeled as \u0026ldquo;increase\u0026rdquo; or \u0026ldquo;decrease,\u0026rdquo; the embedding vector was rotated\u0026thinsp;+\u0026thinsp;90\u0026deg; or \u0026minus;\u0026thinsp;90\u0026deg;, respectively, using Rodrigues\u0026rsquo; rotation formula. This encoding captured the directional semantics of the observed biological changes.\u003c/p\u003e\n\u003cp\u003eGraph construction was performed using PyTorch Geometric (ver. 2.6.1) \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e, in which nodes were defined based on the \u0026ldquo;model main\u0026rdquo; categories (model nodes) and individual \u0026ldquo;target\u0026rdquo; entries (target nodes) extracted from each publication. Within each study, directed edges were established from the model node to each associated target node. Additionally, undirected edges were formed between all target nodes that were co-analyzed within the same study, resulting in a fully connected subgraph among targets. Each edge was encoded as a 1,536-dimensional tensor by concatenating the 768-dimensional feature embeddings of the source and target nodes. To support node-level learning, node embeddings were also aggregated into mean vectors representing the average features across identical model or target nodes, thereby capturing generalizable node representations across the dataset. Full implementation details, including preprocessing scripts and embedding procedures, are available via the GitHub repository (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/fumikawano-lab/Bioreaction-Variation-Network\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eGNN model learning\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eGNN training was conducted using a five-layer architecture (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e). Full implementation details, including preprocessing scripts and model training code, are available via the GitHub repository (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/fumikawano-lab/Bioreaction-Variation-Network\u003c/span\u003e\u003c/span\u003e). To capture Model-to-Target interactions, the first layer employed a multi-head Graph Attention Convolution (GATConv) mechanism \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e to compute attention weights (attn_W) for each target node. Model and target features, initially encoded as 768-dimensional BioBERT-based embeddings, were linearly transformed into 2,048-dimensional representations composed of 8 attention heads with 256 hidden dimensions per head. Three types of attention weights were computed:\u003c/p\u003e\n\u003cul\u003e\n\u003cli\u003e\n\u003cp\u003eattn_m_t: derived from model features and the difference between model and connected target features (i.e., edge features),\u003c/p\u003e\n\u003c/li\u003e\n\u003cli\u003e\n\u003cp\u003eattn_m_n: similarly based on model-to-target differences,\u003c/p\u003e\n\u003c/li\u003e\n\u003cli\u003e\n\u003cp\u003eattn_t: derived from target features and edge features.\u003c/p\u003e\n\u003c/li\u003e\n\u003c/ul\u003e\n\u003cp\u003eThese weights were used to update target node features via a weighted sum:\u003c/p\u003e\n\u003cp\u003etarget_x_new\u0026thinsp;=\u0026thinsp;attn_m_t \u0026times; edge_feature\u0026thinsp;+\u0026thinsp;attn_m_n \u0026times; target_feature\u0026thinsp;+\u0026thinsp;attn_t \u0026times; target_feature\u003c/p\u003e\n\u003cp\u003eSubsequent layers (Layers 2\u0026ndash;5) implemented a multi-layer perceptron (MLP) to capture Target-to-Target interactions. Each layer comprised fully connected (FC) transformations, batch normalization (BN), and rectified linear unit (ReLU) activation.\u003c/p\u003e\n\u003cul\u003e\n\u003cli\u003e\n\u003cp\u003eThe second layer standardized all target node features while preserving dimensionality.\u003c/p\u003e\n\u003c/li\u003e\n\u003cli\u003e\n\u003cp\u003eThe third layer reduced the feature dimensionality from 2,048 to 256 through three sequential FC transformations.\u003c/p\u003e\n\u003c/li\u003e\n\u003cli\u003e\n\u003cp\u003eThe fourth layer was designed to model local topological structure via message passing. Message passing was defined as the difference between a given target node\u0026rsquo;s feature and the average feature of its input neighbors, as derived from the original training data. The mean squared error (MSE) between this reference message passing and the output of the fourth layer was computed as the message passing loss, a component of the total training loss. This reference message passing was then multiplied with each target node\u0026rsquo;s output feature to weight nodes according to their connectivity and propagation potential.\u003c/p\u003e\n\u003c/li\u003e\n\u003cli\u003e\n\u003cp\u003eThe fifth layer modeled pairwise dominance among target nodes using a domination weight (dom_W). The domination weight was computed as the Euclidean distance between the input feature from a specific source node and the mean feature vector of all input features converging on the same target node, including the input from the source node itself. This distance was then scaled by a sigmoid function. A smaller distance indicates that the input feature closely resembles the aggregate input context of the target node, suggesting a stronger upstream (dominant) influence. In the fifth layer, the updated target feature was obtained through two sequential FC\u0026ndash;BN\u0026ndash;ReLU transformations, and trained to approximate the reference feature weighted by dom_W. The MSE between the updated features and the reference was used as domination loss, another component of the total training loss.\u003c/p\u003e\n\u003c/li\u003e\n\u003c/ul\u003e\n\u003cp\u003eModel optimization was performed using Adam optimizer \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e with a learning rate of 0.001, and backpropagation was carried out over 50 or 200 total epochs. To mitigate gradient bias, total loss was computed as a weighted sum of the two loss components: message passing loss \u0026times; 2.0 and domination loss \u0026times; 0.1. Final trained models are publicly available at Google Cloud Storage (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://storage.googleapis.com/skeletal_muscle/sm_v1/gnn_model/gnn_model_final.pt\u003c/span\u003e\u003c/span\u003e for 200-epoch model and gnn_model_final_50epoch.pt for 50-epoch model).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eNetwork inference\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo infer individual-specific hidden pathways, network traversal was performed based on user-provided input data comprising experimental model metadata and observed target parameters. These inputs, formatted consistently with the training data structure, were embedded into 768-dimensional feature vectors using BioBERT, following the same preprocessing pipeline used during model training. To initiate the inference process, cosine similarity was computed between the input model vector and all model embeddings in the training dataset. The top five most similar model nodes were selected, and their directly connected target nodes were identified as primary targets, serving as start nodes for subsequent traversal (Step 1 in Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e). The aim was to reconstruct viable paths from these start nodes to goal nodes, which corresponded to the target features present in the user-provided input. In Step 2, all possible directed paths between the identified start and goal nodes were enumerated based on the GNN-learned edge connectivity between target nodes. This yielded a candidate path space representing potential mechanistic routes through the network. In Step 3, path refinement was performed by selecting the most contextually relevant edges at each intermediate node. When multiple edges existed between the same node pair (i.e., same source and target nodes but differing feature contexts), the edge whose source feature was most similar to the preceding target node, as determined by cosine similarity, was selected. This ensured feature continuity across the reconstructed path. In Step 4, path quality was assessed using a message passing loss. Initial features at the primary nodes were weighted by the learned attention values (attn_W). For each subsequent node, the predicted feature was updated by multiplying the incoming feature by a sigmoid-scaled Euclidean distance derived from the preceding node, thereby propagating context-specific signal from source to goal. The discrepancy between predicted and input-derived goal node features was computed as the message passing loss for each path.\u003c/p\u003e\n\u003cp\u003eSteps 3 and 4 were embedded in a genetic algorithm framework using the DEAP library \u003csup\u003e\u003cspan class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e. Each candidate path generated in Step 2 was treated as an individual. Fitness was evaluated as a weighted combination of message passing loss and diversity, the latter reflecting the uniqueness of nodes and path topology. Path selection employed the Elitist Non-dominated Sorting Genetic Algorithm (NSGA-II), with mutation applied to 20% of individuals per generation. Mutation involved random substitution of an edge within a path, followed by structural repair using the same edge selection strategy described in Step 2 and 3. If the mutation resulted in higher loss, the original path configuration was retained. Given that the initial path candidates were exhaustively constructed and optimized during Steps 2 and 3, the evolutionary loop was iterated twice, solely to confirm that no further improvement in fitness metrics could be achieved. All source code and execution details for the network inference procedure are available at the associated GitHub repository (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/fumikawano-lab/Bioreaction-Variation-Network\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eAnalysis of individualized networks\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo identify individual-specific pathways, inferred GNN outputs for each exercised mouse (n\u0026thinsp;=\u0026thinsp;3) were compared against all non-exercised controls (n\u0026thinsp;=\u0026thinsp;3). For each exercised mouse, individual networks were constructed by aggregating all predicted paths obtained from pairwise comparisons with each control mouse. A common network was then defined as the set of paths shared across all three individual networks, while unique networks were obtained by subtracting the common network from each individual's network. To assess overall reconstruction accuracy, message passing loss was calculated for each goal node and averaged across the network. In addition, among non-primary and non-goal nodes, the most frequently occurring source node across all inferred networks was identified. For this node, the contribution of each connected target node was assessed by referencing the message passing loss of the specific path (i.e., evolutionary algorithm-derived individual) in which the corresponding edge was included. This provided a direct measure of each target node\u0026rsquo;s involvement in the reconstruction of input-derived goal node features within individualized networks.\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cem\u003eOverview of training data\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe training graph was constructed from a total of 27,155 model nodes and 84,723 target nodes. These nodes were interconnected by 383,225 model-to-target edges and 2,475,502 target-to-target edges. Node frequency distributions are summarized in the Supplementary data.xlsx online. Briefly, the most frequently represented experimental models included exercise, high-fat diet feeding, sarcopenia, aging, muscle injury, electrical stimulation, type 2 diabetes, Duchenne muscular dystrophy, diabetes, and isometric contraction. Frequently analyzed target nodes encompassed IL-6, myogenin, MyoD, TNFα, PGC-1α, creatine kinase, insulin, myosin heavy chain, MuRF1, and atrogin1.\u003c/p\u003e\u003cp\u003eDuring model training, the domination loss decreased more rapidly than the message passing loss (see Supplementary Fig. \u003cspan refid=\"MOESM3\" class=\"InternalRef\"\u003eS3\u003c/span\u003e online), a trend that was accounted for by applying differential weighting to these two components in the total loss function. When equal weights were assigned, the message passing loss plateaued early in the training process, limiting further optimization. Although the domination loss approached a minimum around 50 epochs, the message passing loss continued to decline steadily, reducing the difference between the two losses from 0.269 at epoch 50 to 0.228 at epoch 200. These results indicate that training for 200 epochs did not result in overfitting and that both 50-epoch and 200-epoch models were retained for downstream inference to assess potential differences in generalization.\u003c/p\u003e\u003cp\u003e\u003cem\u003eValidation of model with virtual input data\u003c/em\u003e\u003c/p\u003e\u003cp\u003eTo evaluate the inference performance of the trained GNN model, two types of virtual input datasets were constructed that differed only in the biological sex of the subjects, with all other experimental conditions held constant. In the networks inferred from the model trained over 50 epochs, key nodes such as myosin heavy chain I, myonuclei, and VO₂ max were frequently identified in both male and female datasets; however, the surrounding edge structures differed noticeably between sexes (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). Similar nodes were also retrieved from the model trained over 200 epochs, but no clear improvement was observed in the mean distance (average message passing loss) compared to the 50-epoch model, suggesting that inference performance had already stabilized by that point.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cem\u003eNetwork inference with real experimental data\u003c/em\u003e\u003c/p\u003e\u003cp\u003eTo investigate individualized transcriptional responses to exercise, we reanalyzed skeletal muscle RNA sequencing data from exercised and non-exercised mice (n\u0026thinsp;=\u0026thinsp;3 each) previously reported in our study \u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. Differentially expressed genes between exercised and non-exercised conditions were used to construct input features, which were then inferred through the trained GNN models (trained for either 50 or 200 epochs). Because the contextual input for the experimental model (e.g., exercise protocol, species, tissue) was identical across the three exercised mice, the same primary nodes were consistently identified. However, the inferred networks showed marked interindividual variation in the structure and composition of intermediate and goal nodes (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). These differences were reflected in the mean message passing loss, which varied across individuals. Notably, inference using the 200-epoch trained model yielded lower average loss values in all individuals compared to the 50-epoch model, indicating improved network reconstruction performance.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eGiven the superior performance of the 200-epoch model, we proceeded to identify common and unique network paths among individuals using this model. A path was defined as a complete sequence of connected nodes from a primary node to a goal node. A total of 27 paths were commonly found across all three individuals, involving key nodes such as AKT, creatine kinase, FOXO3A, mitochondrial morphology, OPA1, and UQCRC2 as intermediate nodes (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA, B). Unique networks were obtained by subtracting this common network from each individual's network (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC). For example, Exer #1 exhibited 65 unique paths, of which 21 overlapped with Exer #2 and 22 with Exer #3 (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA). Exer #2 and Exer #3 shared 35 paths (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA), indicating that, among the three exercised mice analyzed, these two mice shared similar characteristics in the individual-specific regulatory mechanisms governing transcriptional responses to exercise.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eTo further explore key regulatory factors potentially explaining interindividual variation, we focused on intermediate source nodes, excluding primary and goal nodes, and identified those with the highest frequency across all individuals. UQCRC2 emerged as the most recurrent source node in the inferred networks. This high recurrence indicated that UQCRC2 contributed most consistently to distinguishing the hidden network structures across individuals. Figure\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e illustrates the target nodes connected to UQCRC2 in each individual. UQCRC2 formed unique edges to citrate synthase, RPS6, cytochrome c, ANXA2, UQCRC2, CD36, MYH7, MYHC2A, and SIRT3 exclusively in Exer #1. While Exer #2 and Exer #3 shared several downstream nodes, edges from UQCRC2 to triglycerides and free fatty acids were specific to Exer #2, whereas P70S6K, CKMT2, ACTN3, and MYOZ1 were strongly associated with Exer #3. All edges identified from both common and unique paths are listed in the Supplementary data.xlsx online, along with the message passing loss of the paths in which each edge is included.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eThis study presents a GNN model capable of inferring latent networks that reflect individual-specific responses to physiological stimuli. The model was trained on a dataset constructed from published studies retrieved using the keyword “skeletal muscle”. As of July 2025, approximately 258,000 articles are indexed under this keyword, of which 104,000 are available as free full-text articles via PubMed Central. Our training dataset included 65,096 of these articles, covering more than one-fourth of the accessible literature on skeletal muscle.\u003c/p\u003e\u003cp\u003eThe predominant experimental models extracted from these publications—such as exercise, high-fat diet, sarcopenia, aging, muscle injury, electrical stimulation, type 2 diabetes, and Duchenne muscular dystrophy—indicate that the dataset reflects biologically relevant contexts closely tied to skeletal muscle research. Accordingly, the extracted parameters and target nodes largely comprised skeletal muscle-related molecules and physiological outcomes. Notably, the model was designed to integrate both molecular and physiological parameters. Among the physiological parameters, non-muscle-related parameters such as VO\u003csub\u003e2\u003c/sub\u003e max, blood glucose, and plasma hormone levels were also included, allowing the reconstruction of networks that reflect both intracellular signaling and systemic physiological outcomes. Rather than focusing solely on inferring molecular cascades, the model learned connections between parameters that frequently co-occur across diverse experimental contexts. This architecture enables the exploration of novel, individualized routes from experimental models to observed phenotypes. This capability is likely supported by the model’s global network structure, which connects parameters across distinct publications through shared nodes and feature similarities. As demonstrated by the inference with virtual input data (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e), even a single biological difference, such as sex, led to the identification of markedly distinct hidden networks. These results highlight the model’s ability to adaptively explore individualized paths.\u003c/p\u003e\u003cp\u003eThe model’s inference behavior was influenced by the domination weight learned during training, as edge selection was strongly biased by the learned domination-based relationships between nodes. One key difference between the 50- and 200-epoch models appears to lie in the edge directionality and path structure, which became increasingly refined with extended training. Validation using virtual input data demonstrated that the 50-epoch model was sufficient to reconstruct the network, likely due to the simplicity of the input, which included only two goal nodes (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). In contrast, inference using real biological input data showed improved performance with the 200-epoch model, presumably because the input was more complex, containing unchanged or noisy parameters, thereby requiring deeper learning to refine edge relevance and directionality (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). These observations collectively suggest that the present GNN functions as a scientifically interpretable, small-scale inference model tailored to skeletal muscle biology, with sufficient capacity to support individualized network inference across a tissue-specific corpus of published studies.\u003c/p\u003e\u003cp\u003eInference using real biological input data derived from the tibialis anterior muscles of mice subjected to a single bout of running following four weeks of training identified Uqcrc2 as a key factor potentially explaining the interindividual variation in transcriptional responses to acute exercise (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). UQCRC2 (ubiquinol-cytochrome c reductase core protein 2) is a component of mitochondrial respiratory chain complex III, essential for the formation of mitochondrial supercomplexes involving complexes I, III, and IV \u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e,\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. Missense mutations in UQCRC2 have been linked to impaired mitochondrial respiration and inherited human disorders \u003csup\u003e\u003cspan additionalcitationids=\"CR18\" citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e–\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. Notably, a previous study reported that high volume of high-intensity interval training significantly increased UQCRC2 protein expression in skeletal muscle, contributing to enhanced ATP production through more efficient formation of respiratory complexes \u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. In the individualized networks inferred by the model, UQCRC2 exhibited variable interactions with distinct downstream nodes across mice, suggesting that its local network context contributed differently to each individual’s gene expression response (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). Importantly, the model's use of domination-based weighting does not imply that UQCRC2 lies upstream in a canonical signaling pathway. Rather, the model hypothesizes directional relationships in which UQCRC2 functions as a dominant influence relative to its connected nodes within the inferred context. For instance, CKMT2, ACTN3, and MYOZ1 were uniquely associated with UQCRC2 in the Exer #3 (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e), implying that these factors might be directly or indirectly modulated by UQCRC2, even though the precise molecular mechanisms linking them remain unknown. This interpretive framework reflects the core capability of the present GNN model: to propose biologically plausible, yet hypothetical, directional associations based on learned patterns across multiple studies. While the inferred relationships should be interpreted as hypothesis-generating rather than confirmatory, the model successfully produced skeletal muscle-specific network structures that align with experimental outcomes and capture individual-level diversity. Taken together, these findings support the utility of this GNN as a domain-focused inference tool with strong contextual validity for skeletal muscle biology.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e"},{"header":"Limitations","content":"\u003cp\u003eThis study introduces a GNN model designed to infer bioreaction-variation networks within a skeletal muscle–specific corpus. The primary aim was not to propose a universally optimized model, but rather to demonstrate a customizable model architecture that can be implemented and adapted at a laboratory scale to address specific biological questions. The model presented here serves as a proof-of-concept example of such a design. However, the current implementation has several limitations as shown below\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003eThe model requires input data that represent differential states between experimental and control conditions in order to infer individualized networks. Consequently, datasets comprising only baseline or resting conditions, without any comparative group, are not applicable for inference with this model. Therefore, the extraction of individual-specific features through inference reflects relative interindividual variation within the cohort, including the control subjects used for comparison.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eExperimental models or analytical parameters for which relevant edges are scarce in the training dataset may lack sufficient connectivity to support robust inference. In such cases, the model may fail to generate appropriate or biologically meaningful network predictions.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003eBecause numerical embeddings generated by BioBERT do not directly reflect the magnitude of biological changes in a quantitatively interpretable form, it was necessary to incorporate information on the direction of change, specifically by labeling each relation as either an increase or a decrease. This approach allowed the bioreaction-variation described in the contexts to be translated into vector space. However, capturing such directional changes at the individual level cannot rely on statistical significance testing, which is generally not applicable to single-subject data. In the present study, RNA sequencing data were used, and a twofold change threshold, commonly employed in transcriptomic analyses, was applied to assign directional labels. For other data types or analytical contexts, effect size metrics such as Cohen’s \u003cem\u003ed\u003c/em\u003e may provide a more appropriate basis for estimating individual-level deviations from the comparison group.\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003ch2\u003eCompeting interests\u003c/h2\u003e\u003cp\u003eThe author declares no competing interests.\u003c/p\u003e\u003c/p\u003e\u003cp\u003e\u003ch2\u003eEthics declarations\u003c/h2\u003e\u003cp\u003eNot applicable.\u003c/p\u003e\u003c/p\u003e\u003ch2\u003eFunding\u003c/h2\u003e\u003cp\u003eThe Author received NO FUNDING for this work.\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eFK conceptualized the study, designed the framework of the graph neural network model, wrote all code, validated the model, analyzed the data, drafted the manuscript, and finalized it.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe datasets generated and/or analysed during the current study are available in the GitHub repository, https://github.com/fumikawano-lab/Bioreaction-Variation-Network. The location of each PyTorch graph dataset on private Google Cloud Storage is also provided in the repository.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eBelsky, D. W. et al. Quantification of the pace of biological aging in humans through a blood test, the DunedinPoAm DNA methylation algorithm. \u003cem\u003eeLife\u003c/em\u003e 9, (2020). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.7554/eLife.54870\u003c/span\u003e\u003cspan address=\"10.7554/eLife.54870\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBelsky, D. W. et al. DunedinPACE, a DNA methylation biomarker of the pace of aging. \u003cem\u003eeLife\u003c/em\u003e 11, (2022). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.7554/eLife.73420\u003c/span\u003e\u003cspan address=\"10.7554/eLife.73420\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBamman, M. M., Petrella, J. K., Kim, J. S., Mayhew, D. L. \u0026amp; Cross, J. M. Cluster analysis tests the importance of myogenic gene expression during myofiber hypertrophy in humans. \u003cem\u003eJ. Appl. Physiol.\u003c/em\u003e \u003cb\u003e102\u003c/b\u003e, 2232\u0026ndash;2239. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1152/japplphysiol.00024.2007\u003c/span\u003e\u003cspan address=\"10.1152/japplphysiol.00024.2007\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2007).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eRoss, R. et al. Precision exercise medicine: understanding exercise response variability. \u003cem\u003eBr. J. Sports Med.\u003c/em\u003e \u003cb\u003e53\u003c/b\u003e, 1141\u0026ndash;1153. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1136/bjsports-2018-100328\u003c/span\u003e\u003cspan address=\"10.1136/bjsports-2018-100328\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBonafiglia, J. T. et al. Examining interindividual differences in select muscle and whole-body adaptations to continuous endurance training. \u003cem\u003eExp. Physiol.\u003c/em\u003e \u003cb\u003e106\u003c/b\u003e, 2168\u0026ndash;2176. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1113/EP089421\u003c/span\u003e\u003cspan address=\"10.1113/EP089421\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eXi, X. et al. A mechanism-informed deep neural network enables prioritization of regulators that drive cell state transitions. \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cb\u003e16\u003c/b\u003e, 1284. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41467-025-56475-9\u003c/span\u003e\u003cspan address=\"10.1038/s41467-025-56475-9\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2025).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhang, Y., Feng, X., Wang, Y. \u0026amp; Shi, K. Deep learning powered single-cell clustering framework with enhanced accuracy and stability. \u003cem\u003eSci. Rep.\u003c/em\u003e \u003cb\u003e15\u003c/b\u003e, 4107. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41598-025-87672-7\u003c/span\u003e\u003cspan address=\"10.1038/s41598-025-87672-7\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2025).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eOhsawa, I. \u0026amp; Kawano, F. Chronic exercise training activates histone turnover in mouse skeletal muscle fibers. \u003cem\u003eFASEB journal: official publication Federation Am. Soc. Experimental Biology\u003c/em\u003e. \u003cb\u003e35\u003c/b\u003e, e21453. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1096/fj.202002027RR\u003c/span\u003e\u003cspan address=\"10.1096/fj.202002027RR\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLee, J. et al. BioBERT: a pre-trained biomedical language representation model for biomedical text mining. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cb\u003e36\u003c/b\u003e, 1234\u0026ndash;1240. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1093/bioinformatics/btz682\u003c/span\u003e\u003cspan address=\"10.1093/bioinformatics/btz682\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eFey, M. \u0026amp; Lenssen, J. E. Fast graph representation learning with PyTorch Geometric. \u003cem\u003eICLR\u003c/em\u003e \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.48550/arXiv.1903.02428\u003c/span\u003e\u003cspan address=\"10.48550/arXiv.1903.02428\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2019).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eVeličković, P. et al. \u003cem\u003eGraph Atten. Networks ICLR\u003c/em\u003e, doi:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.48550/arXiv.1710.10903\u003c/span\u003e\u003cspan address=\"10.48550/arXiv.1710.10903\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eZhou, J. et al. Graph neural networks: A review of methods and applications. \u003cem\u003eAI Open.\u003c/em\u003e \u003cb\u003e1\u003c/b\u003e, 57\u0026ndash;81. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.aiopen.2021.01.001\u003c/span\u003e\u003cspan address=\"10.1016/j.aiopen.2021.01.001\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2020).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eKingma, D. P., Ba, J. L. \u0026amp; Adam A method for stochastic optimization. \u003cem\u003eICLR\u003c/em\u003e \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.48550/arXiv.1412.6980\u003c/span\u003e\u003cspan address=\"10.48550/arXiv.1412.6980\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2015).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eDe Rainville, F. M., Fortin, F. A., Gardner, M. A., Parizeau, M. \u0026amp; Gagn\u0026eacute;, C. DEAP: A Python framework for evolutionary algorithms. \u003cem\u003eProceedings of the 14th Annual Conference on Genetic and Evolutionary Computation (GECCO\u003c/em\u003e), 85\u0026ndash;92, ), 85\u0026ndash;92, (2012). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1145/2330784.2330799\u003c/span\u003e\u003cspan address=\"10.1145/2330784.2330799\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2012).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eLetts, J. A., Fiedorczuk, K. \u0026amp; Sazanov, L. A. The architecture of respiratory supercomplexes. \u003cem\u003eNature\u003c/em\u003e \u003cb\u003e537\u003c/b\u003e, 644\u0026ndash;648. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/nature19774\u003c/span\u003e\u003cspan address=\"10.1038/nature19774\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eWu, M., Gu, J., Guo, R., Huang, Y. \u0026amp; Yang, M. Structure of Mammalian Respiratory Supercomplex I(1)III(2)IV(1). \u003cem\u003eCell\u003c/em\u003e \u003cb\u003e167\u003c/b\u003e (e1510), 1598\u0026ndash;1609. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.cell.2016.11.012\u003c/span\u003e\u003cspan address=\"10.1016/j.cell.2016.11.012\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2016).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eBurska, D. et al. Homozygous missense mutation in UQCRC2 associated with severe encephalomyopathy, mitochondrial complex III assembly defect and activation of mitochondrial protein quality control. \u003cem\u003eBiochim. et Biophys. acta Mol. basis disease\u003c/em\u003e. \u003cb\u003e1867\u003c/b\u003e, 166147. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.bbadis.2021.166147\u003c/span\u003e\u003cspan address=\"10.1016/j.bbadis.2021.166147\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGaignard, P. et al. UQCRC2 mutation in a patient with mitochondrial complex III deficiency causing recurrent liver failure, lactic acidosis and hypoglycemia. \u003cem\u003eJ. Hum. Genet.\u003c/em\u003e \u003cb\u003e62\u003c/b\u003e, 729\u0026ndash;731. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/jhg.2017.22\u003c/span\u003e\u003cspan address=\"10.1038/jhg.2017.22\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2017).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eMiyake, N. et al. Mitochondrial complex III deficiency caused by a homozygous UQCRC2 mutation presenting with neonatal-onset recurrent metabolic decompensation. \u003cem\u003eHum. Mutat.\u003c/em\u003e \u003cb\u003e34\u003c/b\u003e, 446\u0026ndash;452. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/humu.22257\u003c/span\u003e\u003cspan address=\"10.1002/humu.22257\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2013).\u003c/span\u003e\u003c/li\u003e\u003cli\u003e\u003cspan\u003eGranata, C. et al. High-intensity training induces non-stoichiometric changes in the mitochondrial proteome of human skeletal muscle without reorganisation of respiratory chain content. \u003cem\u003eNat. Commun.\u003c/em\u003e \u003cb\u003e12\u003c/b\u003e, 7056. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41467-021-27153-3\u003c/span\u003e\u003cspan address=\"10.1038/s41467-021-27153-3\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2021).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"interindividual variation, artificial intelligence, machine learning, deep learning","lastPublishedDoi":"10.21203/rs.3.rs-7239284/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7239284/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eInterindividual variation in biological responses to physiological stimuli is a widely recognized phenomenon. However, effective computational tools for identifying the individual-specific mechanisms remain limited. We present a graph neural network (GNN) model designed to infer hidden molecular and physiological relationships underlying such variation in experimental biological data. To ensure applicability at a laboratory scale, the model was trained on a domain-specific corpus constructed from approximately 65K published studies containing the keyword \u0026ldquo;skeletal muscle\u0026rdquo;. The architecture comprises five layers with a multi-head attention mechanism and a multi-layer perceptron, enabling the model to capture both local topological features and directional dominance between connected nodes. The GNN was trained to learn relationships from experimental models to target features, as well as among target features. Using real experimental input consisting of differential gene expression data from mouse skeletal muscle subjected to acute exercise, the model successfully inferred individualized networks, identifying both common and unique paths across individuals based on input experimental context. These results demonstrate the model\u0026rsquo;s capacity to extract interpretable, individual-specific biological connectivity patterns. The proposed framework serves as a proof of concept for customizable, context-based GNN inference designed to address biological variation at the individual level.\u003c/p\u003e","manuscriptTitle":"A graph neural network model for inferring interindividual variation from experimental biological data","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-08-07 08:57:21","doi":"10.21203/rs.3.rs-7239284/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-08-13T14:50:39+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-08-04T17:12:35+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"105393831155331358518030284464731401696","date":"2025-08-04T16:25:46+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-08-04T09:22:49+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"108444051692051054489144750257343265936","date":"2025-08-04T07:15:43+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"300973840835648595773995261156302981170","date":"2025-08-04T06:43:50+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-08-04T06:38:07+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-08-04T06:17:04+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2025-07-31T17:44:28+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-07-30T12:35:44+00:00","index":"","fulltext":""},{"type":"submitted","content":"Scientific Reports","date":"2025-07-30T12:32:15+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"7141c406-e2b3-4414-85fb-1ac88c3137d4","owner":[],"postedDate":"August 7th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":52738540,"name":"Biological sciences/Computational biology and bioinformatics"},{"id":52738541,"name":"Physical sciences/Mathematics and computing"},{"id":52738542,"name":"Biological sciences/Neuroscience"},{"id":52738543,"name":"Biological sciences/Systems biology"}],"tags":[],"updatedAt":"2025-11-17T16:00:55+00:00","versionOfRecord":{"articleIdentity":"rs-7239284","link":"https://doi.org/10.1038/s41598-025-23320-4","journal":{"identity":"scientific-reports","isVorOnly":false,"title":"Scientific Reports"},"publishedOn":"2025-11-12 15:57:03","publishedOnDateReadable":"November 12th, 2025"},"versionCreatedAt":"2025-08-07 08:57:21","video":"","vorDoi":"10.1038/s41598-025-23320-4","vorDoiUrl":"https://doi.org/10.1038/s41598-025-23320-4","workflowStages":[]},"version":"v1","identity":"rs-7239284","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7239284","identity":"rs-7239284","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}