Mechanochemical modeling of exercise

This is an uncorrected proof. Figures Abstract Skeletal muscle displays remarkable plasticity, adapting its size and strength in response to mechanical loading, especially, from exercise. This process, known as hypertrophy, is fundamental to athletic training and rehabilitation, but is challenging to quantitatively predict due to its multifactorial, multiscale nature. Specifically, skeletal muscle hypertrophy results from an integration of macroscopic mechanical stimuli with the intracellular signaling pathways that govern muscle growth. In this work, we present a multiscale computational model that mechanistically integrates these mechanical and biochemical stimuli and offers a framework for predicting the outcomes of different types of exercise on skeletal muscle growth. The framework couples a transversely isotropic hyperelastic model for tissue-level mechanics with a system of ordinary differential equations representing the IGF1-AKT-mTOR-FOXO signaling pathway, a key regulator of protein synthesis and degradation. We link these scales using a volumetric growth model, where the signaling dynamics inform a growth tensor that drives changes in muscle cross-sectional area. This approach enables the simulation of long-term muscle adaptation, providing a mechanistic tool to investigate how different exercise protocols lead to macroscopic hypertrophy. Simulations from our model capture the temporal dynamics of hypertrophy under varying load protocols and highlight how feedback between protein synthesis and muscle growth regulates the dose-response relationship to prevent unbounded growth. Using muscle geometries derived from the Visible Human dataset, we study how human variations in muscle geometry affect hypertrophy. Finally, we demonstrate that the mechanochemical coupling between muscle geometry and signaling not only predicts macroscopic shape changes but also provides buffering from local signaling heterogeneity. Ultimately, this framework offers a predictive computational tool for optimizing training regimens and understanding the multiscale determinants of muscle adaptations. Author summary Skeletal muscle naturally adapts its size and strength in response to physical activity, a process that is essential for both athletic training and rehabilitation. However, predicting how a specific training routine will alter a muscle’s shape remains difficult. This is because muscle growth relies on a complex, multiscale chain of events, where physical forces from movements trigger chemical signals that slowly rebuild the tissue over time. To better understand this process, we developed a computational model that links these underlying cellular signals directly to the macroscopic physical changes in the muscle. When simulating different exercise routines in realistic human muscle geometries, we found that the complex internal fiber architecture of the muscle causes it to grow unevenly. We also observed that the physical continuity of the tissue acts as a mechanical buffer, helping to coordinate noisy, localized signaling into smooth overall growth. Ultimately, this framework provides a foundational first step towards understanding how and where muscle adds mass, laying the groundwork for more advanced predictive tools. Citation: Devold IS, Rognes ME, Rangamani P (2026) Mechanochemical modeling of exercise-induced skeletal muscle hypertrophy. PLoS Comput Biol 22(8): e1014691. https://doi.org/10.1371/journal.pcbi.1014691 Editor: Emma Lejeune, Boston University, UNITED STATES OF AMERICA Received: February 14, 2026; Accepted: August 10, 2026; Published: August 24, 2026 Copyright: © 2026 Devold et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: Data and source code are available via Zenodo: Devold IS. Supplementary Code: Mechanochemical Modeling of Exercise-Induced Skeletal Muscle Hypertrophy; 2026. Zenodo. https://doi.org/10.5281/zenodo.17902802. Funding: This work was supported in part by the Wu Tsai Human Performance Alliance at the University of California, San Diego (to P.R.). I.S.D is supported by the Simula–UCSD–University of Oslo Research and PhD training (SUURPh) program, an international collaboration in computational biology and medicine funded by the Norwegian Ministry of Education and Research. M.E.R. acknowledges support from Stiftelsen Kristian Gerhard Jebsen via the K.G. Jebsen Centre for Brain Fluid Research, Wellcome via Award 313298/Z/24/Z (Next-generation simulation and learning in imaging-based biomedicine), and from the Research Council of Norway (RCN). Simulations were performed on the Experimental Infrastructure for Exploration of Exascale Computing (eX3), which is financially supported by the RCN under contract 270053. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: I have read the journal’s policy and the authors of this manuscript have the following competing interests: P.R. is a consultant for Simula Research Laboratory in Oslo, Norway and receives income. The terms of this arrangement have been reviewed and approved by the University of California, San Diego in accordance with its conflict-of-interest policies. Introduction Skeletal muscle is remarkably plastic, constantly adapting its size and structure to mechanical stimuli resulting from different movements. This plasticity is not only critical for muscle’s ability to respond to the physical demands placed on it [1,2], but also critical to athletic training, rehabilitation, and overall health. Different types of exercise, such as high-intensity sprints, high-load resistance training, or long-duration endurance training, each induce distinct patterns of muscle activation that differ in duration, amplitude, and frequency (Fig 1A). Although the effects of exercise on skeletal muscle function are empirically well-documented [3], the underlying connections between molecular mechanisms and mechanical responses driving these adaptations remain poorly understood [4]. This is largely because they represent a profoundly multiscale phenomenon that requires an integration of tissue-level mechanics and subcellular biochemical signal transduction that has proven challenging to model. (A) Different exercise types are characterized by distinct patterns of muscle activation, leading to physiological adaptations and ultimately affecting muscle performance. Here, “progress” refers to improvements in physical performance (strength), which is driven initially by rapid neural adaptations and later by structural muscle growth. (B) In this work, we focus on the IGF1-AKT signaling pathway, which is a key regulator of muscle hypertrophy. Mechanical stimuli through exercise activate IGF1, which in turn activates AKT. AKT promotes protein synthesis via mTOR and inhibits protein degradation via FOXO. (C) Schematic of the proposed multiscale model, coupling a signaling model with tissue-level mechanics to predict exercise-induced muscle growth. Abbreviations: CSA: cross-sectional area; IGF1: insulin-like growth factor 1; AKT: protein kinase B; mTOR: mechanistic target of rapamycin; FOXO: forkhead box transcription factor. Figure partially created in BioRender. Rangamani, P. (2026) [https://BioRender.com/wfzzhb7]. The multiscale nature of muscle spans both spatial and temporal scales [4]. Spatially, muscle tissue is organized hierarchically, from the whole organ down to fascicles, individual muscle fibers, myofibrils, and finally to the sarcomere – the fundamental 2–3 µm contractile unit composed of actin and myosin filaments [5,6]. Temporally, muscle contraction is driven by rapid electrochemical events. This includes millisecond-timescale calcium transients triggering the ATP-fueled cross-bridge cycle, which in turn generates force and causes the relative sliding of the sarcomere filaments [7]. Over longer timescales, that span weeks to months, the cumulative effects of these contractions trigger intricate signaling pathways that regulate protein synthesis and degradation [8,9] and changes to metabolic pathways [10] ultimately remodeling the tissue structure. Mechanistically linking an exercise protocol to a long-term adaptation requires integrating our understanding of both muscle mechanics and cellular physiology. From a mechanical perspective, foundational models of cellular-level force generation [11,12] have, with modern advances in medical imaging and computational resources, evolved into detailed three-dimensional models [13]. These models typically employ a hyperelastic material framework to accurately capture the large, non-linear deformations that muscle tissue undergoes during activity [14–16]. More recently, some models have begun to bridge scales by coupling cell-level models [17] with tissue-level mechanics [18–20]. From a physiological perspective, there are multiple signaling pathways at play in the muscle [8]. Among these, the IGF1-AKT-mTOR-FOXO signaling pathway has been identified as a central regulator of load-induced skeletal muscle hypertrophy [21–23]. In this cascade, mechanical loading through exercise stimulates the release of insulin-like growth factor (IGF1), which activates the protein kinase B (AKT) (Fig 1B). AKT, in turn, displays a dual response: it promotes protein synthesis via mechanistic target of rapamycin (mTOR) and inhibits protein degradation by suppressing the forkhead box transcription factor (FOXO), suggesting the presence of an incoherent feedforward loop [24,25]. The balance between these mTOR and FOXO-driven mechanisms dictates the net change in muscle protein and, consequently, muscle growth. Computational models shed critical insights into how exercise can regulate muscle growth at different length and time scales. For short time-scale events, systems biology approaches have extensively mapped the calcium dynamics underlying excitation-contraction coupling [17,26], with recent work clarifying the role of store-operated Ca2+ entry (SOCE) in force modulation [7], the link between mechanical loading and metabolic sensors such as AMP-activated protein kinase (AMPK) [27], and the bidirectional coupling between exercise and mitochondrial dynamics [28]. However, such signaling models typically represent the muscle tissue as a homogeneous compartment, predicting hypertrophy as a scalar quantity without accounting for the complex spatial distribution of stress and strain within the 3D tissue. Conversely, tissue-level mechanical models excel at resolving these local deformations using hyperelastic formulations. However, they predominantly operate under the assumption of mechanical equilibrium with fixed material reference states. By treating the tissue as a static material, these models neglect the dynamic, biochemical adaptation of the tissue constituents that occurs in response to mechanical loading over time. Bridging these domains requires overcoming a fundamental timescale mismatch: instantaneous mechanical stimuli trigger signaling events that accumulate over weeks to drive morphological change, which in turn alters the mechanical environment. The mathematical theory of volumetric growth [29–31] offers a rigorous kinematic framework to resolve this multiscale challenge. By decomposing the deformation gradient into elastic and growth contributions, this theory allows biologically driven mass accumulation from tissue growth to be coupled directly to continuum mechanics. The volumetric growth framework has previously found applications in several biological tissues [32,33], including cardiac [34,35] and brain tissue [36]. In the context of skeletal muscle, Villota-Narvaez et al. [37] found that a multiscale mechanobiological model, which integrates the IGF1-AKT signaling pathway with continuum mechanics via a growth tensor, accurately predicts muscle hypertrophy (CSA) results under various training conditions. However, their framework focused on idealized geometries and did not fully implement the cumulative growth kinematics required to robustly track residual stresses and shape evolution over long-term training protocols. In this paper, we seek to address the following challenges: How can we simulate growth kinematics over timescales spanning several weeks to months while properly accounting for evolving residual stresses? How can we incorporate the effect of complex, non-uniform myofiber architectures found in realistic muscle anatomies and their influence on local growth? And finally, how can we leverage advances in computational mechanics to generate a predictive framework for muscle responses to different exercise schemes? To this end, we present a multiscale computational framework that couples tissue-level mechanics with the IGF1-AKT signaling pathway (Fig 1C). We ground our approach in the mathematical theory of volumetric growth [29,31], which allows biologically driven mass accumulation to be coupled directly to continuum mechanics. While our work builds on the mechanobiological model of Villota-Narvaez et al. [37], it extends the computational framework in two important directions. First, we adopt the cumulative growth framework of Goriely and Ben Amar [38], which ensures that each growth increment is applied to a stress-free configuration, providing a robust method for handling long-term adaptation. Second, we apply this coupled model to anatomically realistic muscle geometries derived from the Visible Human dataset [39–41], demonstrating the versatility of the approach. Using this framework, we demonstrate that high-intensity, high-frequency exercise protocols drive the most pronounced hypertrophy, reveal that complex muscle architectures result in non-uniform growth patterns, and that introducing heterogeneity in the signaling dynamics disrupts the overall growth. This integrated approach demonstrates how protocol-dependent signaling drives 3D tissue adaptation in realistic geometries, highlighting its predictive potential for both rehabilitation and peak performance sports. Methods Stress-strain relationships Muscle tissue undergoes large deformations and exhibits a nonlinear stress-strain relationship under physiological loading. Its bundle-like structure results in a transversely isotropic material response, with distinct mechanical properties along and across the fiber directions [42]. Accordingly, we model the mechanical behavior of muscle using a transversely isotropic hyperelastic framework [19]. In this framework, the stress response is derived from a strain energy density function, W, which depends on the tissue deformation. To define this relationship, we first introduce key kinematic quantities. Let denote a fixed reference domain (open, bounded) with coordinates and boundary . We let denote the displacement from this reference configuration to the current configuration. The deformation gradient is given by , and the corresponding right Cauchy-Green deformation tensor is . The local volume change is measured by the Jacobian . The strain energy function for a transversely isotropic material is defined through the invariants of : the first invariant and the fourth pseudo-invariant . Finally, the fiber stretch is given by . The stress experienced by the material is quantified with respect to the reference configuration by the first and second Piola–Kirchhoff stress tensors, and , which are derived from a strain energy function as and , respectively, and are related through the deformation gradient as . In the constitutive model presented by Röhrle et al. [19], the second Piola–Kirchhoff stress tensor is given by the sum of isotropic, anisotropic and active contributions: (1)In (1), the isotropic term defines a Mooney–Rivlin material with material constants c1 and c2 (Table 1) and hydrostatic pressure p. The term accounts for passive fiber anisotropy, while represents active fiber contraction. The force-length relationships, fpassive and factive, depend on the fiber stretch and follow a standard normalized model [43] where the maximal active force is generated at an optimal fiber stretch, (Fig 2A) [14,18]. Specifically, we define the passive and active components by the following piecewise exponential and polynomial functions, respectively [14,18]: (2a)(2b)where the parameters and are given in Table 1, while and are chosen to ensure continuity and differentiability. Upon activation (), the active stress component generates contraction along the fiber direction (Fig 2B, 2F). (A) The passive and active fiber forces are dependent on stretch . (B) A muscle contraction, showing the inactive () and active () states. (C) Illustration of one week of the three exercise protocols used in growth simulations. (D) The coupling scheme between signaling and growth mechanics used in our computational framework. The exercise model determines an incremental growth tensor , which updates the mechanical model through the cumulative growth framework. The resulting cross-sectional area feeds back to modulate the signaling model via the protein synthesis rate . (E) Schematic of the cumulative growth framework, based on the multiplicative decomposition of deformation, . Between growth steps, computational unloading (, green pathway) is required to avoid applying new deformation to a residually stressed state (red arrow). (F) Example geometries showing the computational mesh, surfaces for boundary conditions (), and prescribed fiber architecture. Mechanical equilibrium equations To determine the muscle tissue’s mechanical state for a given activation level , we solve for the displacement-pressure pair satisfying the governing equations of static equilibrium defined relative to the reference domain : (3a)(3b)(3c)(3d)Here, (3a) represents mechanical equilibrium in the absence of body forces, (3b) enforces incompressibility, and is the outward unit normal to the boundary . The Robin boundary includes the two ends of the muscle, denoted as and , respectively. On these specific surfaces, we apply the nominal traction via the Robin boundary condition (3c). We assume that the physical Cauchy traction , which is defined relative to the deformed boundary surface with unit normal , is negatively proportional to the displacement in the normal direction with proportionality constant : This condition models the muscle connecting to elastic springs with stiffness k (Fig 2B), allowing contraction while stabilizing the nonlinear solver. We then define as the pull-back of to the reference configuration: where we have used Nanson’s formula to map between area in the reference and deformed configurations. The remaining boundaries are traction-free (3d). Signaling model The IGF1-AKT-mTOR-FOXO signaling pathway is a central regulator of muscle mass, balancing protein synthesis and degradation [8]. To model the biochemical drivers of hypertrophy, we adopt the signaling model from Villota-Narvaez et al [44], which in turn builds upon a simplified representation of the pathway by Schiaffino and Mammucari [22] (Fig 1B). The model consists of a system of ordinary differential equations (ODEs) describing the pathway’s key components. Let I, A, F, and M denote the normalized populations of IGF1, AKT, FOXO, and mTOR, respectively. Their dynamics are modeled using a Lotka–Volterra framework that captures the activating and inhibiting interactions between them through intrinsic growth rates (), self-inhibition rates (), and coupling strengths () (Table 2): (4a)(4b)(4c)(4d)The rate of change of the myofibril population N reflects a balance between mTOR-driven synthesis and FOXO-driven degradation. We first define the unconstrained signaling rate, f(F, M), which captures these contributions when the signals exceed their respective thresholds, and : (5)To ensure the myofibril population remains within an admissible range , we further impose: (6)This balance is driven by the protein synthesis rate (Fig C in S1 Appendix) and degradation rate , with contributions effective only when mTOR and FOXO exceed their respective thresholds, and . Exercise is modeled as a stimulus that modulates the intrinsic growth rates of both IGF1 and AKT. First, the IGF1 growth rate, , is elevated during exercise according to the relation (7)where is a scaling factor, is the exercise intensity factor, is the muscle’s maximum force capacity, and denotes the baseline rate (Table 2). The stimulus is subject to a physiological delay, which we represent by in the governing equation for IGF1 (4a) and set to h. Second, the intrinsic growth rate of AKT, , responds to exercise with a transient increase defined by (8)In this expression, is a delay parameter, a time-scaling factor, and marks the start time of the most recent exercise session. This formulation creates a transient increase in AKT’s growth rate that decays over several hours, while the max operation ensures does not fall below its baseline level, . Exercise protocols We simulated three distinct exercise regimens to test protocol-dependent adaptation (Fig 2C). These were: a high-intensity, high-frequency protocol (1h at 80% intensity, 3x/week, MWF) based on DeFreitas et al. [45]; a high-intensity, low-frequency protocol (1h at 80% intensity, 1x/week); and a low-intensity, medium-frequency protocol (4h at 20% intensity, every three days). The intensity factor was set to 0.8 and 0.2 for the high-intensity and low-intensity regimes, respectively. For all protocols, there is an initial week without exercise to establish baseline conditions, followed by eight weeks of the specified exercise regimen. Initial conditions The initial conditions () for the signaling model represent a homeostatic, pre-exercise steady state. To determine these values, we simulated the ODE system (4) with baseline parameter values and no external stimulus () until a stable equilibrium was reached. These equilibrium values, listed in Table 2, were used for all subsequent simulations. Spatial variations in signaling To capture cell-to-cell variability and represent the muscle as a heterogeneous tissue, we introduce spatial variation in the signaling model. Specifically, spatial heterogeneity and inter-individual variability are emulated by perturbing selected signaling parameters locally. Each baseline parameter p0 (Table 2) is replaced by a spatially varying parameter where is a random field sampled from a uniform distribution on (e.g., for ±10% variation). Perturbations were drawn independently per mesh cell to represent within-muscle heterogeneity, and sampled once per simulation to simulate between-subject variability. Sensitivity analysis We performed a Sobol sensitivity analysis [46] using the SALib Python library [47,48] to identify the most influential parameters in the signaling model. The final myofibril population (N) after a 15-day simulated exercise protocol was selected as the primary output metric, serving as a computationally efficient proxy for muscle growth. For the analysis, each parameter was sampled from a uniform distribution spanning 50% to 200% of its baseline value. To ensure a robust assessment, the analysis included first- and second-order effects with 10,000 samples, resulting in 300,000 total model evaluations. Model coupling To capture the interplay between molecular signaling and tissue mechanics during skeletal muscle hypertrophy, we couple the subcellular signaling pathway with the macroscopic mechanical model. It is important to note that in the current framework, growth is driven directly by the signaling output of the ODE model rather than by local mechanical stimuli such as computed tissue stress or strain. Consequently, aggregate phenomenological parameters governing the whole-muscle exercise stimulus, such as the maximum force capacity , are decoupled from the maximum isometric stress parameter that strictly dictates tissue-level continuum mechanics. We implement this specific form of bidirectional coupling through two primary mechanisms: (1) the signaling dynamics inform a volumetric growth tensor that drives changes in muscle geometry, and (2) the resulting macroscopic change in muscle size provides feedback to the signaling pathway by modulating a key parameter – the protein synthesis rate (Fig 2D). These two coupling mechanisms are detailed in the following sections. From signaling to macroscopic growth through the growth tensor Molecular signals are translated into macroscopic growth by linking the myofibril population N(t) from the signaling model to the tissue’s deformation. This coupling is formalized using the theory of volumetric growth, where the total deformation gradient is decomposed into a growth part and an elastic part , such that [29,31]. In this decomposition, the growth tensor represents the local (potentially incompatible) addition of mass driven by biology, while is the subsequent elastic deformation required to maintain mechanical equilibrium and ensure the tissue remains a compatible, continuous body. Importantly, the biological tissue is assumed incompressible under elastic loading (3b). Therefore, the incompressibility constraint is imposed strictly on the elastic tensor, . This allows the total tissue volume to increase dynamically in response to exercise-driven growth. To determine , we first calculate a scalar growth multiplier, G(t). The core assumption linking N(t) to the continuum growth is that the myofibril population is directly proportional to the muscle cross-sectional area (CSA), , via a constant , such that . The explicit formulation for calculating in the deformed geometry is developed later in the text (13). From the change in myofibrils predicted by the signaling model over a time step , the growth multiplier is defined as: (9)where is the myofibril population at the previous time step. This scalar multiplier then defines the anisotropic growth tensor : (10)This formulation ensures that growth occurs exclusively in the plane perpendicular to the fiber direction , consistent with an increase in CSA [32]. Crucially, the 3D continuum framework translates the zero-dimensional scalar output of the signaling model into physical, non-uniform tissue deformations determined by local fiber architecture. Feedback from mechanics to signaling The bidirectional coupling is completed by a feedback mechanism where the mechanical state of the muscle modulates the signaling pathway. Concretely, the protein synthesis rate is set to depend on the current CSA of the muscle . We define the normalized CSA as , where is the initial area. The synthesis rate is then updated according to a Hill-type inhibitory function: (11)With a Hill coefficient of n = 4 and a half-max constant of , this function establishes homeostatic control. The synthesis rate returns to its baseline value when , decreases during hypertrophy () to limit excessive growth, and increases during atrophy () to stimulate growth back towards the baseline. This bounded, Hill-type regulation provides a simple, yet physiologically plausible feedback loop for muscle size adaptation (Figs D,E in S1 Appendix). Implementation of cumulative growth To model the long-term effects of incrementally applied growth, we adopt the cumulative growth framework of Goriely and Ben Amar [38]. A key requirement of this framework is that each new growth increment, , must be applied to a stress-free configuration. Because the tissue is generally under residual stress from previous growth steps, this requires computationally “unloading” the accumulated elastic deformation before applying the next increment (Fig 2E). After k growth steps, the total deformation is decomposed into a total elastic tensor and a total growth tensor , such that . These total tensors are updated recursively at each step k: (12a)(12b)(12c)In this formulation, the incremental growth tensor is determined by the signaling model via (10). equation (12b) integrates this new growth by first computationally unloading the previous state. Specifically, the inverse of the accumulated elastic tensor, , maps the tissue from its current, residually-stressed configuration back to a virtual stress-free state. The new growth increment is then applied to this relaxed state, and the kinematics are pushed forward again via . Finally, the finite element solver computes mechanical equilibrium for this updated total growth tensor , yielding the new incremental elastic tensor , which updates the total elastic deformation via (12c). A step-by-step pseudo-code of this numerical routine is provided in Algorithm A in S1 Appendix. Muscle geometries Mesh generation. Simulations were performed on three distinct types of computational geometries: a cylindrical geometry, an idealized fusiform muscle, and a set of anatomically realistic muscle models (Fig 2F). The idealized fusiform geometry was created to represent a simplified, yet physiologically relevant, muscle shape. The tetrahedral mesh was generated using Gmsh [49] and consists of 1997 cells. The realistic muscle geometries were derived from segmentations of the Visible Human Dataset [39,40]. From these surface representations, tetrahedral meshes were generated using fTetWild [50], resulting in cell counts between 6707 and 9668. Muscle fiber architecture. Defining the muscle fiber architecture, which sets the local anisotropy , is essential for all geometry types. For the cylinder, fibers simply run parallel along the long axis. In the idealized fusiform model, fibers were defined analytically via a spline curve. For the realistic geometries, fiber directions were estimated by solving a Stokes flow problem [51]. This involved setting the muscle origin and insertion as inflow and outflow boundaries, applying slip conditions on the remaining surface, and normalizing the resulting velocity field to define the local fiber vectors (Fig 2F). Muscle cross-sectional area. The cross-sectional area (CSA) was computed by integrating over a transverse reference surface intersecting the muscle at the midpoint of its longitudinal axis. Nanson’s formula [52] provides the mapping from the reference to the deformed area element for this integration: (13)For the idealized geometry, was defined as an internal circular surface located at the geometric midpoint and explicitly included during mesh generation in Gmsh. For each realistic geometry, a unique reference cross-section was defined using principal component analysis (PCA): the primary principal component determined the longitudinal muscle axis, and was taken as the plane passing through the geometric centroid and normal to this axis. The displacement field was interpolated onto this plane to compute the deformed area. Numerical solution schemes and verification To simulate the coupled mechanics-signaling system, we use a quasi-static approach, justified by the slow timescale of hypertrophy compared to that of elastic deformation [38]. The simulation advances through a series of exercise and rest periods, using a time step of h. Each period begins by solving the signaling ODE system to determine the rate of change in the myofibril population. This growth is then applied incrementally in a nested loop, where each growth step consists of constructing a new growth tensor, solving for mechanical equilibrium, and finally, closing the feedback loop by updating the protein synthesis rate based on the current cross-sectional area (Algorithm A in S1 Appendix). The computational model was implemented in Python, and the simulation code is openly available [53]. The signaling ODE system was integrated using the ‘RK45’ method from SciPy’s solve_ivp function [54], with a local time step of dt = 0.05h. For the spatially varying signaling models, however, a custom implementation of the fourth-order Runge–Kutta (RK4) method was used. We use the finite element method to numerically solve the mechanical equilibriumnonlinear equations (3), using the FEniCSx finite element software [55]. Relative to each of the computational meshes, we define the finite element space as the space of continuous piecewise quadratic 3-vector fields and Q as the space of continuous piecewise linears. We then introduce the test functions and . At each iteration, we solve the following system of nonlinear equations: find the approximate displacement and pressure solving (14a)(14b)where is the first Piola–Kirchhoff stress tensor, with defined in (1). The nonlinear system (14) was solved using a Newton solver. To verify the implementation, we performed a series of numerical tests, including evaluating the correctness of the numerical solution algorithms for the spatially varying ODE systems (Fig A in S1 Appendix) and a mesh and time step convergence test for the full growth simulation (Fig B in S1 Appendix). Results Signaling model links exercise protocol to distinct dynamic responses Before implementing the full multiscale framework, we first characterized the dynamics of the core signaling model in isolation. We simulated the model’s response to three distinct exercise protocols (Fig 2C): (1) a high-intensity, high-frequency protocol (MWF); (2) a high-intensity, low-frequency protocol (weekly); and (3) a low-intensity, medium-frequency protocol (every three days). The model exhibited a dynamically consistent response across all protocols (Fig 3A-3E). Each exercise bout triggered a rapid, transient spike in the upstream signals IGF1 and AKT, with a peak observed approximately 12 hours post-stimulus (Fig 3A-3B). The amplitude of these spikes was directly proportional to the exercise intensity. The upstream activation drove the expected opposing response in the downstream effectors: mTOR levels (promoting synthesis) increased, while FOXO levels (promoting degradation) were suppressed (Fig 3C-3D). Importantly, the system operates on two distinct timescales. IGF1 and AKT act as fast responders, returning to baseline between each exercise period. In contrast, FOXO and mTOR display a wave summation effect, wherein the responses to repeated stimuli accumulate and oscillate around an apparent steady state. The magnitude of this shift is protocol-dependent and directly determines the final myofibril density N (Fig 3E). (A–E) Time courses of the normalized populations (A) IGF1, (B) AKT, (C) FOXO, (D) mTOR, and (E) myofibrils N for three different exercise protocols. Insets in (A) and (B) display a zoomed-in 24-hour time window to highlight the rapid transient spikes in IGF1 and AKT signaling. (F) Sobol sensitivity analysis quantifying the influence of each model parameter on the final myofibril number obtained using a standardized daily exercise protocol. The plot shows the mean Sobol index versus its standard deviation for each parameter. Marker colors identify the associated signaling molecule, with split colors indicating coupling parameters. To identify which model parameters were most critical in controlling this response, we performed a Sobol sensitivity analysis. The analysis revealed that the final myofibril population was most sensitive to parameters governing the downstream effector, FOXO (Fig 3F). Specifically, the intrinsic growth () and self-inhibition () rates of FOXO were the most influential, followed by the corresponding rates for mTOR ( and ) and AKT ( and ). This sensitivity arises because the core signaling network functions as an incoherent feedforward loop [24,25]: the upstream signal (AKT) simultaneously promotes a synthesis pathway (mTOR) and suppresses a degradation pathway (FOXO). Consequently, the qualitative outcome – whether the muscle undergoes hypertrophy or atrophy – is not determined by the exercise protocol alone but depends critically on the kinetic balance between these opposing branches. The high sensitivity to FOXO parameters indicates that and effectively set a threshold for growth. Small perturbations to these rates can abruptly toggle the net protein turnover from positive to negative, even under identical exercise stimuli. Comparing the net growth response between the regimens, the high-intensity, high-frequency (MWF) protocol induced the largest increase in the myofibril population (N) due to the summation effects described above (Fig 3E). However, all exercise protocols predicted net growth relative to the sedentary baseline, confirming that even lower-frequency stimuli provide a sufficient perturbation to shift the net protein balance. With these protocol-dependent trajectories of N established, we next used them as the driving input for the full 3D framework to investigate how biochemically driven growth manifests as mechanical remodeling across idealized muscle geometries. Coupled model predicts protocol-dependent hypertrophy in an idealized muscle geometry While biochemical signaling and tissue mechanics are often modeled in isolation, physiological growth is an inherently coupled phenomenon. Accordingly, we next integrated the signaling model with the 3D hyperelastic framework, initially simulating long-term adaptation in an idealized fusiform muscle geometry (Fig 4A). (A) Fiber orientation a0 computed from a spline expression. (B) Mechanical state variables (displacement, fiber stretch , and fiber stress) at zero and full activation (). (C–E) Predicted time courses of cross-sectional area, muscle volume, and the protein synthesis rate in response to three distinct exercise protocols. (F) Final muscle configurations after the simulation period for each protocol, with the growth-induced displacement amplified 10x for visualization. First, we characterized the muscle’s baseline mechanical behavior in the absence of volumetric growth. A simulated contraction from zero to full activation ( from 0 to 1) produced heterogeneous distributions of displacement, stress and fiber stretch along the muscle geometry (Fig 4B). Muscle ends were allowed to partially contract towards the midplane. The displacement magnitudes were largest at the ends of the muscle, consistent with this shortening behavior. Tensile fiber-direction Cauchy stress peaked at the surface of the muscle center, transitioning to negative (compressive) values internally, whereas the fiber stretch was relatively uniformly distributed, exhibiting only slightly higher values in the central region. Having established that the mechanical framework gives the expected response for an idealized muscle geometry, we next combined the signaling model with the muscle growth model. We simulated muscle adaptation over a nine-week period for the three different exercise protocols. The coupled model successfully translated the protocol-dependent signaling dynamics into different macroscopic growth in an exercise-dependent manner. Consistent with the myofibril-level predictions, the high-intensity, high-frequency MWF protocol induced the most significant increase in both cross-sectional area (CSA) and total muscle volume, followed by the every-three-days and weekly protocols (Fig 4C-4E). Furthermore, we also identified the critical role of homeostatic control generated by the feedback between signaling and mechanics. As the muscle CSA increased in response to training, the mechanical feedback loop progressively downregulated the protein synthesis rate, , from its baseline value (Fig 4C, 4E). The strength of this feedback is dependent on the exercise modality – higher exercise frequency leads to a sharper drop in the protein synthesis rate, naturally limiting the rate of further hypertrophy (Fig 4D). The macroscopic effect of this adaptation is visualized in Fig 4F, which shows the final deformed muscle shapes for each protocol (growth amplified tenfold for visualization). This result is a direct consequence of the anisotropic growth tensor included in the model (10), which mechanistically links the signaling output to mass addition in the plane perpendicular to the fiber axis. Overall, these simulations suggest that while any exercise stimulus promotes growth relative to the baseline, the mechanical feedback loop prevents a linear dose-response, ensuring the higher-frequency training leads to self-limiting adaptation rather than unbounded hypertrophy. Anatomically realistic models show protocol-dependent adaptation and non-uniform growth Thus far, we used idealized geometries to investigate how muscle growth responds to different exercise modalities. To determine how complex, subject-specific anatomy influences this response, we next applied the coupled framework to a cohort of 12 anatomically realistic muscle geometries derived from the Visible Human Dataset [39,40]. The cohort consisted of bilateral (left and right) meshes for both a male and female subject, yielding four unique geometries for each of the three selected muscle types: the biceps femoris long head (BFLH), semitendinosus (ST), and tibialis anterior (TA) (Fig 5A). Functionally, the BFLH and ST act as hamstring knee flexors and hip extensors, while the TA serves as a primary ankle dorsiflexor. These muscles are heavily recruited during common activities such as running, squatting, and lunging, making them highly relevant for investigating exercise-induced adaptation. Prior to simulation, we characterized the baseline geometries of the cohort, confirming that male muscles displayed larger initial volumes and cross-sectional areas compared to their female counterparts (Fig 5B). Furthermore, to capture the internal architecture essential for anisotropic growth, we generated subject-specific fiber arrangements for each mesh using Stokes flow simulations (see Methods, Fig 5C, and Fig F in S1 Appendix). (A) Anatomical placement of the semitendinosus (ST), biceps femoris long head (BFLH), and tibialis anterior (TA) muscles, visualized within the Visible Human Female. (B) Comparison of computed muscle volumes and cross-sectional areas for the twelve meshes derived from the bilateral ST, BFLH, and TA muscles of the Visible Human Female and Male. (C) Fiber arrangements generated using Stokes flow simulations within meshes derived from the surface geometries shown in (A). Labels indicate gender (F: Female, M: Male) and laterality (L: Left, R: Right). We first examined the model’s predictions for a representative case: the female left biceps femoris long head (BFLH). The 3D visualizations of the muscle following the nine-week training protocols reveal a distinct, protocol-dependent hypertrophic response (Fig 6A). Consistent with the high-intensity nature of the stimulus, the MWF protocol produced the most substantial visible thickening of the muscle belly, followed by the every-three-days and weekly protocols. Our simulations show that the complex, subject-specific architecture results in irregular tissue expansion, in contrast to the symmetric radial scaling observed in the idealized model. While the idealized geometry distributes growth evenly around the central axis (Fig 4F), the realistic geometry exhibits spatially heterogeneous hypertrophy, with peak displacement localized to specific regions of the muscle belly. To understand the regulation underlying this growth, we tracked the temporal evolution of the protein synthesis rate, (Fig 6B). The simulations show a progressive downregulation of from its baseline value as the tissue adapts. This decay reflects the homeostatic feedback loop: as the muscle grows, the mechanical stimulus for further synthesis is reduced, naturally limiting the adaptation rate over time. (A) Final muscle configurations of the female left biceps femoris long head (BFLH) following the simulation period for each exercise protocol. Growth-induced displacements are amplified 10x for visualization. (B) Predicted temporal evolution of the protein synthesis rate () in response to three distinct loading protocols (shown for the female left BFLH geometry). (C–D) Predicted final CSA and volume after three protocols in the BFLH, ST, and TA muscles. Error bars indicate standard error. Expanding this analysis to the full cohort of 12 geometries, we found that these protocol-dependent trends were robust across all subjects and muscle types. For every geometry tested, the high-frequency, high-intensity MWF protocol induced the most pronounced hypertrophy, followed by the every-three-days and weekly protocols (Fig 6C-6D). A key finding from the cohort analysis was that the relative increase in total muscle volume was consistently larger than the relative increase in cross-sectional area (CSA) measured at the muscle’s midplane. For instance, the MWF protocol induced an average volume increase of approximately 8–9% (Fig 6D), while the corresponding CSA increase was only 3–6% (Fig 6C). This discrepancy arises because the complex architecture drives non-uniform growth, where mass is added asymmetrically along the muscle. Consequently, site-specific CSA measurements, a common metric in clinical and experimental studies [56], may systematically underestimate the extent of tissue adaptation. The variability was notably more pronounced in the normalized CSA measurements (Fig 6C) compared to the more tightly clustered volume measurements (Fig 6D). This reinforces that CSA is a localized metric, highly sensitive to anatomical placement and local geometric irregularities, whereas volume more robustly captures the global adaptation. Within the scope of our simulations, we did not see any muscle type specific or sex specific differences in response to exercise. In conclusion, these results demonstrate that tissue geometry and fiber architecture play a fundamental role in the spatial distribution of hypertrophy. Because the model defines growth as mass addition transverse to the local fiber axis, the subject-specific fiber orientation directly dictates the direction of expansion at every point. Consequently, complex fiber arrangements result in non-uniform shape changes that cannot be fully captured by idealized models, highlighting the need for anatomical fidelity in predictive modeling of muscle adaptation. Tissue mechanics integrates heterogeneous signaling into smooth deformation Finally, to mimic cell-to-cell variability in muscle tissue, we investigated how inherent spatial variation in biochemical signaling affects the mechanochemical response to exercise. Moving beyond the assumption of homogeneous tissue properties, we conducted an ensemble analysis on the female left BFLH geometry using the high-intensity MWF protocol. This specific combination was selected as the “maximal response” case to determine if robust hypertrophy persists even under noisy signaling conditions. In these simulations, all rate parameters associated with IGF1, AKT, FOXO, and mTOR were simultaneously perturbed by ±1% using a uniform random distribution at each mesh cell. Crucially, given the incoherent feed-forward loop structure of the signaling network, such parameter fluctuations have the potential to qualitatively alter the local outcome, flipping the kinetic balance from net synthesis to degradation. The time courses of the local ODE states reveal the model’s sensitivity to cellular variability (Fig 7A). While the ensemble’s mean signaling response remained consistent with the uniform baseline, the variance was significantly amplified in the downstream regulators, FOXO and mTOR. This amplification reflects a key biological feature: because downstream molecules integrate signals over time, they accumulate and magnify small upstream fluctuations. Consequently, even under the optimal high-intensity (MWF) stimulus, local outcomes vary significantly. Although the average response is positive, specific regions may experience blunted growth or even atrophy due to unfavorable parameter combinations. (A) Predicted time courses of the local ODE states (IGF1, AKT, FOXO, mTOR, and myofibril population N). Heterogeneity was introduced by sampling rate parameters (, , ) from uniform distributions within of their baseline values. Gray lines and shading indicate the mean standard deviation across five samples; the dashed blue line represents the baseline uniform case. (B) Predicted evolution of normalized cross-sectional area (CSA), normalized volume, and the protein synthesis rate under the MWF protocol (colors as in A). (C) Spatial maps of final myofibril density (N), local volume change (J), and growth-induced displacement magnitude (amplified 10x) for the uniform case versus a representative heterogeneous sample. At the macroscopic level, the uniform baseline case consistently achieved more total growth than the mean of the heterogeneous ensemble, as measured by both normalized volume and CSA (Fig 7B). This indicates that the introduction of signaling noise creates a net drag on the system, where the cumulative effect of suboptimal local responses slightly dampens the aggregate adaptation compared to the idealized, uniform condition. This difference in growth is reflected in the protein synthesis rate, . Since the uniform case achieves a larger increase in CSA, it triggers a correspondingly larger decrease in compared to the heterogeneous average. Furthermore, this analysis confirmed our previous observation that variability was consistently greater for normalized CSA than for total volume. Because CSA is a localized slice, it is highly susceptible to specific regional fluctuations, whereas total volume integrates over the entire domain, effectively smoothing out the local noise to provide a more robust metric of adaptation. The 3D spatial maps reveal a striking contrast between the biological drive and the mechanical response (Fig 7C). The heterogeneity in the myofibril population (N) is closely reflected in the speckled, fragmented map of local growth (J). However, the resulting displacement magnitude shows a smoother and more coordinated pattern. This phenomenon demonstrates a mechanical buffering, where the global equilibrium constraints of the tissue force the continuum to deform smoothly. Effectively, the elastic matrix averages out the highly variable local expansion rates, preventing individual cells from expanding or shrinking in isolation. However, this buffering comes at a cost. The mismatch between the jagged local growth drive (J) and the smooth deformation gradient implies the generation of internal residual stresses. It is possible that regions of muscle tissue attempting to grow rapidly are physically constrained by their slower-growing neighbors, creating a mechanical drag that not only dampens the aggregate expansion but may also generate localized shear concentrations prone to micro-damage. Discussion In this work, our goal was to develop a multiscale computational framework to quantify the effects of different exercise modalities on realistic human muscle geometries. This framework successfully integrated IGF1-AKT signaling dynamics with 3D hyperelastic tissue mechanics to model long-term skeletal muscle hypertrophy. We first established that the readouts of the signaling pathway are highly sensitive to the exercise protocol, with the high-intensity, high-frequency regimen inducing the most significant growth. The architecture of the signaling network resembles an incoherent feedforward loop [24,25]. We identified FOXO and mTOR as the critical points governing the balance between protein synthesis and degradation for muscle hypertrophy, though further experimental validation is required to confirm this. The same stimulus resulting in a balance of synthesis and degradation is a common feature found in many signaling pathways [57–59] and is important for maintaining homeostasis. Applying the coupled model to anatomically realistic human geometries confirmed this protocol-dependent adaptation, highlighting that the feedback between mechanical deformation and signaling is essential for regulating the magnitude of hypertrophy and stabilizing the growth process. Our simulations also revealed that the relative increase in total muscle volume consistently exceeds the increase in cross-sectional area (CSA) in realistic geometries. Finally, our analysis of spatial heterogeneity revealed that while local volume expansion rates are highly localized by biological variability, the tissue-wide integration of mechanics acts as a buffer, forcing the overall displacement and resulting bulk shape change to remain relatively smooth and coordinated. While our current framework lacks the matched longitudinal data required for a rigorous quantitative validation, its outputs fall within physiologically plausible orders of magnitude when compared to anatomical and physiological benchmarks. First, the baseline muscle geometries derived from the Visible Human dataset are consistent with MRI-derived population averages. For instance, Handsfield et al. [60] reported mean volumes of 206.5 cm3 for the BFLH, 186.0 cm3 for the ST, and 135.2 cm3 for the TA. These values compare favorably with our model’s initial volumes of 252.4 cm3, 190.1 cm3, and 132.4 cm3, respectively, supporting the anatomical realism of the meshes but not validating the growth model itself. Regarding the link between mechanical workload and biological response, our model successfully predicted the dose-dependent nature of hypertrophy. We observed that myofibril accretion is highly sensitive to the exercise regimen, with high-frequency, high-intensity protocols inducing the steepest monotonic increase. This prediction aligns with the systematic review by Schoenfeld et al [61], which identifies a significant positive relationship between weekly training frequency (in the range 1-3x/week) and muscle growth. Building on this biochemical foundation, the predicted magnitude of macroscopic hypertrophy falls within a plausible range compared to existing experimental literature. For example, our simulation of the DeFreitas (MWF) protocol on realistic geometries predicted a total volume increase of nearly 9%. Although not a direct validation, as it involves a different muscle group and experimental context, this is on the same order of magnitude as the 10.1-11.3% volume increase reported by Kubo et al [62] in the pectoralis major following 10 weeks of training. Interestingly, the associated increase in mid-belly cross-sectional area was more conservative ( 4–6%) than the range observed in comparable resistance training studies. Specifically, our predictions are slightly lower than the 9.6% increase in thigh muscle CSA reported by DeFreitas et al [45], the 10% increase found by Damas et al. [63] in vastus lateralis after 10 weeks of resistance training, and the 7.4% increase reported by Seynnes et al [64] in the quadriceps femoris following a five-week high-intensity resistance training program. Unlike idealized geometries where volume scales predictably with CSA (Fig 4C, 4D), anatomical muscles exhibit spatially heterogeneous growth, such that local CSA changes may not directly reflect whole-muscle hypertrophy. Thus, our results suggest that single-slice CSA measurements may underrepresent the total hypertrophic response in complex muscle shapes. This may partly explain why our predicted CSA increases are smaller than those reported experimentally. However, the discrepancy may also reflect deficiencies in the current framework. The underprediction of CSA may stem from the highly reduced nature of the signaling model or the use of generalized material parameters that do not capture the exact biomechanical behavior of the specific muscle types studied. A rigorous quantitative validation of the coupled model will require longitudinal datasets suitable for parameter calibration and independent testing, which were not available for the present study. While the proposed framework captures key mechanistic links between signaling and tissue mechanics, we note some limitations and opportunities for future computational development. First, the signaling model is a reduced representation of the underlying biochemical network and omits additional pathways and fast-timescale dynamics. A complex signaling model could capture more biochemical interactions, but at the cost of higher computing demands and the introduction of parameters that are hard to estimate [65]. A future modeling direction could involve coupling detailed signaling models [28,66,67] with multicompartmental, mixed-dimensional PDE growth models [68] to represent different length scales from subcellular organelles to muscle cross-section. Our framework also combines biochemical and mechanical models developed from different sources. Although this integration enables a comprehensive computational framework, establishing parameter provenance and confidence across coupled models remains a long-standing challenge in computational and systems biology that will require concerted effort across the field to resolve [69,70]. Second, we model growth as transverse expansion driven by a global signaling output, rather than by localized mechanical stimuli (e.g., tissue stress or strain). As a result, the aggregate 0D exercise parameters (, ) are completely decoupled from the parameters of the continuum model (, ). Relating these time-averaged, macroscopic variables to local mechanotransduction pathways remains an open challenge for future modeling. Furthermore, our model does not account for along-fiber remodeling, known as sarcomerogenesis [71]. This means the model captures how the tissue thickens (hypertrophy) but misses how fibers lengthen to adapt to chronic stretch, an important aspect of tissue remodeling. Accurate representation of either growth mechanism ultimately depends on the underlying fiber architecture, since local fiber orientation dictates the axis of both force generation and tissue expansion. Here, fiber paths were estimated using Stokes flow, which provides a mathematical idealization that enforces a generalized alignment from origin to insertion but cannot capture more complex architectures such as the bipennate structure of the tibialis anterior (TA). These generated fiber fields should therefore be viewed as computational proxies rather than high-fidelity anatomical reconstructions. Future improvements could incorporate liquid-crystal models, as recently demonstrated for cardiac muscle [72], while the framework is also readily extensible to patient-specific fiber architectures derived from diffusion tensor imaging (DTI) as such data become more available. Nevertheless, the current cohort provides a valuable testbed for demonstrating the framework’s sensitivity to 3D spatial heterogeneity and its potential for personalized simulations. The mechanical model also remains simplified in several aspects. Incorporating muscle-specific material parameters as experimental measurements become available, for instance, in the force-length relationships, could improve predictive capability. Likewise, the boundary conditions employ linear elastic springs primarily to ensure computational stability rather than to faithfully represent tendon mechanics. Future work could instead incorporate patient-specific tendon and aponeurosis geometries, enabling more realistic nonlinear tissue compliance. Finally, assigning independent random values to each element makes the spatial heterogeneity length scale mesh-dependent. We adopt this uncorrelated distribution as a theoretical stress test to demonstrate the capacity for mechanical buffering. Since actual biological heterogeneity likely manifests as smoother, spatially correlated fields, we expect the qualitative buffering effect to remain, although its quantitative characteristics may differ. This study establishes a modular computational infrastructure for linking cellular biology to organ-level morphology. A key feature of this mechanochemical framework is its extensibility. The reduced IGF1-AKT model serves as an adaptable component that can be substituted with more expansive systems biology models to capture phenomena such as metabolic fatigue, calcium handling, or gene regulation. On the mechanical side, the kinematic description of growth can be extended to include sarcomerogenesis [73]. Incorporating this mechanism in future iterations could enable the differentiation between concentric-induced thickening and eccentric-induced fascicle lengthening, an important capability for analyzing diverse training modalities. A critical next step will be closing the loop on mechanochemical coupling by transitioning from global signaling outputs to local mechanical drivers. Moreover, developing the computational strategies necessary to quantify accumulated residual stresses, alongside evaluating additional metrics like pennation angle, curvature, or physiological CSA, could provide a more comprehensive characterization of the structural adaptations than simple macroscopic volume and anatomical CSA. Beyond mathematical extensions, our framework lays the groundwork for integrating the next generation of complex anatomical data. With the capability to model growth on realistic geometries established, there is a need for larger geometric cohorts and subject-specific fiber architectures derived from DTI. Our platform is ready to incorporate such datasets as they become available, allowing for the rigorous characterization of how anatomical variability dictates functional adaptation. Ultimately, by bridging the gap between molecular drivers and macroscopic function, our computational framework provides the foundation for simulating individual-specific scenarios, from optimizing athlete peak performance to designing targeted rehabilitation protocols for muscle atrophy. Supporting information S1 Appendix. Fig B. Growth simulation verification. Fig C. Parameter fit of the protein synthesis rate . Fig D. Comparison of Hill-type, linear, and no feedback. Fig E. Illustration of reduced growth rate due to feedback. Fig F. Illustration of muscle fiber directions. Algorithm A in S1 Appendix. Cumulative growth implementation. https://doi.org/10.1371/journal.pcbi.1014691.s001 (PDF) Acknowledgments The authors are grateful to Henrik Finsberg and Jørgen Dokken for their assistance with the FEniCSx finite element implementation, and to Halvor Herlyng for generously sharing his implementation of Stokes flow with slip boundary conditions. Declaration of conflicts: P.R. is a consultant for Simula Research Laboratory in Oslo, Norway and receives income. The terms of this arrangement have been reviewed and approved by the University of California, San Diego in accordance with its conflict-of-interest policies. References - 1. Egan B, Sharples AP. Molecular responses to acute exercise and their relevance for adaptations in skeletal muscle to exercise training. Physiol Rev. 2023;103(3):2057–170. pmid:36395350 - 2. Roberts MD, McCarthy JJ, Hornberger TA, Phillips SM, Mackey AL, Nader GA, et al. Mechanisms of mechanical overload-induced skeletal muscle hypertrophy: current understanding and future directions. Physiol Rev. 2023;103(4):2679–757. pmid:37382939 - 3. Joyner MJ, Coyle EF. Endurance exercise performance: the physiology of champions. J Physiol. 2008;586(1):35–44. pmid:17901124 - 4. Wisdom KM, Delp SL, Kuhl E. Use it or lose it: Multiscale skeletal muscle adaptation to mechanical stimuli. Biomechanics and Modeling in Mechanobiology. 2015;14(2):195–215. - 5. Rassier DE. Muscle Biophysics. Springer New York. 2010. https://doi.org/10.1007/978-1-4419-6366-6 - 6. Marieb E, Hoehn K. Muscles and muscle tissue. Human anatomy and physiology. 11 ed. Pearson Education. 2019. p. 280–8. - 7. Francis EA, Hamid J, Kumar A, Rangamani P. Systems modeling reveals that store-operated calcium entry modulates force and fatigue during exercise. bioRxiv. 2025. https://doi.org/10.1101/2025.05.22.655415 - 8. Sandri M. Signaling in muscle atrophy and hypertrophy. Physiology (Bethesda). 2008;23:160–70. pmid:18556469 - 9. Blaauw B, Schiaffino S, Reggiani C. Mechanisms modulating skeletal muscle phenotype. Compr Physiol. 2013;3(4):1645–87. pmid:24265241 - 10. Hargreaves M, Spriet LL. Skeletal muscle energy metabolism during exercise. Nat Metab. 2020;2(9):817–28. pmid:32747792 - 11. Hill AV. The heat of shortening and the dynamic constants of muscle. Proceedings of the Royal Society of London Series B - Biological Sciences. 1938;126(843):136–95. - 12. Huxley AF. Muscle Structure and Theories of Contraction. Progress in Biophysics and Biophysical Chemistry. 1957;7:255–318. - 13. Wakeling JM, Febrer-Nafría M, De Groote F. A Review of the Efforts to Develop Muscle and Musculoskeletal Models for Biomechanics in the Last 50 Years. Journal of Biomechanics. 2023;155:111657. - 14. Blemker SS, Pinsky PM, Delp SL. A 3D model of muscle reveals the causes of nonuniform strains in the biceps brachii. J Biomech. 2005;38(4):657–65. pmid:15713285 - 15. He X, Taneja K, Chen J-S, Lee C-H, Hodgson J, Malis V, et al. Multiscale modeling of passive material influences on deformation and force output of skeletal muscles. Int J Numer Method Biomed Eng. 2022;38(4):e3571. pmid:35049153 - 16. Almonacid JA, Domínguez-Rivera SA, Konno RN, Nigam N, Ross SA, Tam C, et al. A Three-Dimensional Model of Skeletal Muscle Tissues. SIAM J Appl Math. 2024;84(3):S538–66. - 17. Shorten PR, O’Callaghan P, Davidson JB, Soboleva TK. A mathematical model of fatigue in skeletal muscle force contraction. J Muscle Res Cell Motil. 2007;28(6):293–313. pmid:18080210 - 18. Röhrle O, Davidson JB, Pullan AJ. Bridging Scales: A Three-Dimensional Electromechanical Finite Element Model of Skeletal Muscle. SIAM J Sci Comput. 2008;30(6):2882–904. - 19. Röhrle O, Davidson JB, Pullan AJ. A physiologically based, multi-scale model of skeletal muscle structure and function. Frontiers in Physiology. 2012;3. - 20. Heidlauf T, Röhrle O. A multiscale chemo-electro-mechanical skeletal muscle model to analyze muscle contraction and force generation for different muscle fiber arrangements. Frontiers in Physiology. 2014;5. - 21. Bodine SC, Stitt TN, Gonzalez M, Kline WO, Stover GL, Bauerlein R, et al. Akt/mTOR pathway is a crucial regulator of skeletal muscle hypertrophy and can prevent muscle atrophy in vivo. Nat Cell Biol. 2001;3(11):1014–9. pmid:11715023 - 22. Schiaffino S, Mammucari C. Regulation of skeletal muscle growth by the IGF1-Akt/PKB pathway: insights from genetic models. Skelet Muscle. 2011;1(1):4. pmid:21798082 - 23. Sandri M, Barberi L, Bijlsma AY, Blaauw B, Dyar KA, Milan G, et al. Signalling pathways regulating muscle mass in ageing skeletal muscle. The role of the IGF1-Akt-mTOR-FoxO pathway. Biogerontology. 2013;14(3):303–23. - 24. Hart Y, Antebi YE, Mayo AE, Friedman N, Alon U. Design principles of cell circuits with paradoxical components. Proc Natl Acad Sci U S A. 2012;109(21):8346–51. pmid:22562798 - 25. Hart Y, Alon U. The utility of paradoxical components in biological circuits. Mol Cell. 2013;49(2):213–21. pmid:23352242 - 26. Senneff S, Lowery MM. Effects of extracellular potassium on calcium handling and force generation in a model of excitation-contraction coupling in skeletal muscle. J Theor Biol. 2021;519:110656. pmid:33667541 - 27. Linden-Santangeli N, Zhang J, Kramer B, Rangamani P. Systems modeling and uncertainty quantification of AMP-activated protein kinase signaling. NPJ Syst Biol Appl. 2025;11(1):113. pmid:41087369 - 28. Khalilimeybodi A, Qiao L, Leung A, McCulloch AD, Schenk S, Rangamani P. Systems modelling of mitochondrial dynamics in different exercise regimes. J Physiol. 2026;604(15):6369–90. pmid:42385159 - 29. Rodriguez EK, Hoger A, McCulloch AD. Stress-dependent finite growth in soft elastic tissues. J Biomech. 1994;27(4):455–67. pmid:8188726 - 30. Lubarda VA, Hoger A. On the mechanics of solids with a growing mass. International Journal of Solids and Structures. 2002;39(18):4627–64. - 31. Goriely A. The mathematics and mechanics of biological growth. Springer New York. 2017. https://doi.org/10.1007/978-0-387-87710-5 - 32. Ambrosi D, Ben Amar M, Cyron CJ, DeSimone A, Goriely A, Humphrey JD, et al. Growth and remodelling of living tissues: perspectives, challenges and opportunities. J R Soc Interface. 2019;16(157):20190233. pmid:31431183 - 33. Kuhl E. Growing matter: a review of growth in living systems. J Mech Behav Biomed Mater. 2014;29:529–43. pmid:24239171 - 34. Göktepe S, Abilez OJ, Kuhl E. A generic approach towards finite growth with examples of athlete’s heart, cardiac dilation, and cardiac wall thickening. Journal of the Mechanics and Physics of Solids. 2010;58(10):1661–80. - 35. Göktepe S, Abilez OJ, Parker KK, Kuhl E. A multiscale model for eccentric and concentric cardiac growth through sarcomerogenesis. J Theor Biol. 2010;265(3):433–42. pmid:20447409 - 36. Budday S, Steinmann P, Kuhl E. The role of mechanics during brain development. J Mech Phys Solids. 2014;72:75–92. pmid:25202162 - 37. Villota-Narvaez Y, Garzón-Alvarado DA, Röhrle O, Ramírez-Martínez AM. Multi-Scale Mechanobiological Model for Skeletal Muscle Hypertrophy. Frontiers in Physiology. 2022;13:899784. - 38. Goriely A, Ben Amar M. On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity. Biomech Model Mechanobiol. 2007;6(5):289–96. pmid:17123061 - 39. Andreassen TE, Hume DR, Hamilton LD, Walker KE, Higinbotham SE, Shelburne KB. Three Dimensional Lower Extremity Musculoskeletal Geometry of the Visible Human Female and Male. Sci Data. 2023;10(1):34. pmid:36653365 - 40. Andreassen TE, Hume DR, Hamilton LD, Higinbotham SE, Shelburne KB. Automated 2D and 3D Finite Element Overclosure Adjustment and Mesh Morphing Using Generalized Regression Neural Networks. Medical Engineering & Physics. 2024;126(1):104136. - 41. Spitzer VM, Whitlock DG. The Visible Human Dataset: the anatomical platform for human simulation. Anat Rec. 1998;253(2):49–57. pmid:9605360 - 42. Morrow DA, Haut Donahue TL, Odegard GM, Kaufman KR. Transversely isotropic tensile material properties of skeletal muscle tissue. J Mech Behav Biomed Mater. 2010;3(1):124–9. pmid:19878911 - 43. Zajac FE. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering. 1989;17(4):359–411. - 44. Villota-Narvaez Y, Garzon-Alvarado DA, Ramirez-Martinez AM. A dynamical system for the IGF1-AKT signaling pathway in skeletal muscle adaptation. Biosystems. 2021;202:104355. pmid:33453318 - 45. DeFreitas JM, Beck TW, Stock MS, Dillon MA, Kasishke PR 2nd. An examination of the time course of training-induced skeletal muscle hypertrophy. Eur J Appl Physiol. 2011;111(11):2785–90. pmid:21409401 - 46. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer Physics Communications. 2010;181(2):259–70. - 47. Herman J, Usher W. SALib: An open-source Python library for Sensitivity Analysis. JOSS. 2017;2(9):97. - 48. Iwanaga T, Usher W, Herman J. Toward SALib 2.0: Advancing the accessibility and interpretability of global sensitivity analyses. Socio-Environmental Systems Modelling. 2022;4:18155. - 49. Geuzaine C, Remacle JF. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering. 2009;79(11):1309–31. - 50. Hu Y, Schneider T, Wang B, Zorin D, Panozzo D. Fast tetrahedral meshing in the wild. ACM Trans Graph. 2020;39(4). - 51. Varvik J, Besier TF, Handsfield GG. Computational fluid dynamics simulations for 3D muscle fiber architecture in finite element analysis: Comparisons between computational fluid dynamics and diffusion tensor imaging. Int J Numer Method Biomed Eng. 2021;37(12):e3521. pmid:34411442 - 52. Holzapfel GA. Nonlinear solid mechanics: A continuum approach for engineering. Sussex: John Wiley & Sons, LTD. 2000. - 53. Devold IS. Supplementary code: mechanochemical modeling of exercise-induced skeletal muscle hypertrophy. Zenodo. 2026. https://doi.org/10.5281/zenodo.17902802 - 54. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17(3):261–72. pmid:32015543 - 55. Baratta IA, Dean JP, Dokken JS, Habera M, Hale JS, Richardson CN. DOLFINx: The next generation FEniCS problem solving environment. Zenodo. 2023. https://doi.org/10.5281/zenodo.10447666 - 56. Lanza MB, Martins-Costa HC, De Souza CC, Lima FV, Diniz RCR, Chagas MH. Journal of Biomechanics. 2022;132:110956. - 57. Qiao L, Getz M, Gross B, Tenner B, Zhang J, Rangamani P. Spatiotemporal orchestration of calcium-cAMP oscillations on AKAP/AC nanodomains is governed by an incoherent feedforward loop. PLoS Comput Biol. 2024;20(10):e1012564. pmid:39480900 - 58. Tenner B, Getz M, Ross B, Ohadi D, Bohrer CH, Greenwald E. Spatially Compartmentalized Phase Regulation of a Ca2 -cAMP-PKA Oscillatory Circuit. eLife. 2020. https://doi.org/10.7554/eLife.55013 - 59. Rangamani P, Levy MG, Khan S, Oster G. Paradoxical signaling regulates structural plasticity in dendritic spines. Proc Natl Acad Sci U S A. 2016;113(36):E5298-307. pmid:27551076 - 60. Handsfield GG, Meyer CH, Hart JM, Abel MF, Blemker SS. Relationships of 35 lower limb muscles to height and body mass quantified using MRI. J Biomech. 2014;47(3):631–8. pmid:24368144 - 61. Schoenfeld BJ, Ogborn D, Krieger JW. Effects of Resistance Training Frequency on Measures of Muscle Hypertrophy: A Systematic Review and Meta-Analysis. Sports Med. 2016;46(11):1689–97. pmid:27102172 - 62. Kubo K, Ikebukuro T, Yata H. Effects of 4, 8, and 12 Repetition Maximum Resistance Training Protocols on Muscle Volume and Strength. J Strength Cond Res. 2021;35(4):879–85. pmid:32304514 - 63. Damas F, Phillips SM, Lixandrão ME, Vechin FC, Libardi CA, Roschel H, et al. Early resistance training-induced increases in muscle cross-sectional area are concomitant with edema-induced muscle swelling. Eur J Appl Physiol. 2016;116(1):49–56. pmid:26280652 - 64. Seynnes OR, de Boer M, Narici MV. Early skeletal muscle hypertrophy and architectural changes in response to high-intensity resistance training. J Appl Physiol (1985). 2007;102(1):368–73. pmid:17053104 - 65. Linden NJ, Kramer B, Rangamani P. Bayesian parameter estimation for dynamical models in systems biology. PLoS Comput Biol. 2022;18(10):e1010651. pmid:36269772 - 66. Fowler A, Knaus KR, Khuu S, Khalilimeybodi A, Schenk S, Ward SR, et al. Network model of skeletal muscle cell signalling predicts differential responses to endurance and resistance exercise training. Exp Physiol. 2024;109(6):939–55. pmid:38643471 - 67. Coccimiglio IF, Clarke DC. ADP is the dominant controller of AMP-activated protein kinase activity dynamics in skeletal muscle during exercise. PLoS Comput Biol. 2020;16(7):e1008079. pmid:32730244 - 68. Francis EA, Laughlin JG, Dokken JS, Finsberg HNT, Lee CT, Rognes ME, et al. Spatial modeling algorithms for reactions and transport in biological cells. Nat Comput Sci. 2025;5(1):76–89. pmid:39702839 - 69. Budde K, Smith J, Wilsdorf P, Haack F, Uhrmacher AM. Relating simulation studies by provenance-Developing a family of Wnt signaling models. PLoS Comput Biol. 2021;17(8):e1009227. pmid:34351901 - 70. Cowan AE, Mendes P, Blinov ML. ModelBricks-modules for reproducible modeling improving model annotation and provenance. NPJ Syst Biol Appl. 2019;5:37. pmid:31602314 - 71. Andrews MH, S AP, Gurchiek RD, Pincheira PA, Chaudhari AS, Hodges PW, et al. Multiscale hamstring muscle adaptations following 9 weeks of eccentric training. J Sport Health Sci. 2025;14:100996. pmid:39461588 - 72. Barnafi NA, Osses A. Modeling of cardiac fibers as oriented liquid crystals. Computer Methods in Applied Mechanics and Engineering. 2025;436:117710. - 73. Zöllner AM, Abilez OJ, Böl M, Kuhl E. Stretching skeletal muscle: chronic muscle lengthening through sarcomerogenesis. PLoS One. 2012;7(10):e45661. pmid:23049683