Introduction

Sarcomeres are the basic contractile cytoskeletal units of striated muscle, including cardiac muscle. In most species, sarcomeres are ∼2 µm in length and consist of highly ordered bipolar bundles of myosin motor proteins interdigitated with actin filaments, tied together in the Z-bands. (Fig. 1A). Force generation by myosin interacting with actin is regulated by intracellular Ca²⁺ and fueled by chemical energy (ATP)¹. Sarcomeres are serially connected into myofibrils that span the length of the cell. A cardiomyocyte (CM) contains tens of parallel myofibrils. Overall CM contractility is a mesoscopic process that emerges from a very large number of coupled, non-linear, stochastic, elementary force generators. Despite extensive knowledge of single-molecule motor kinetics^2,3 and whole-muscle mechanics^4,5, this intermediate mesoscopic scale—where tens of sarcomeres interact within a myofibril—remains largely unexplored in cardiac muscle. Sarcomeres are mechanically coupled in series, yet each is also an active, excitable unit undergoing rapid, transient activation-relaxation cycles. Do these coupled sarcomeres maintain synchronized motion, or do stochastic fluctuations and non-equilibrium forces drive heterogeneous behavior?

Figure 1:

Recent studies, utilizing in vivo imaging and advanced in vitro models, have begun to reveal notable heterogeneity and beat-to-beat variability at the sarcomere level — even as the average myofibril motion remains smooth and regular^6-8. However, these observations have been limited to murine systems (≈5-8 Hz beating rate vs. ∼1 Hz in humans). While sarcomere non-uniformity has a long history in skeletal muscle research^9-12, prevailing models of cardiac contraction have historically assumed that sarcomeres behave as uniform, synchronized units, with force generation determined primarily by the steady-state force-length relationship—an idea rooted in the Frank-Starling mechanism and length-dependent activation^4,13.

Existing attempts to explain cardiac sarcomere heterogeneity remain largely qualitative, invoking static structural differences between sarcomeres—e.g., in strength or resting length—especially pronounced in cardiomyopathy^6-8,14. Such models propose deterministic mechanisms: ‘weak’ sarcomeres are consistently stretched by ‘strong’ neighbors¹⁴, diastolic length variations dictate subsequent systolic interactions and potential role-switching between sarcomeres^6,8, or inherent length variability enables recruitment of more contracting sarcomeres upon cell lengthening through length-dependent activation⁷.

Crucially, these deterministic, primarily length-centric views neglect velocity-dependent active force generation, friction effects, and non-linearities. By assuming a continuous balance between passive and active forces dictated by steady-state relationships, they ignore the non-equilibrium nature of the beating heart, where adaptation kinetics and viscous drag are critical. Consequently, unraveling these dynamics requires novel approaches that resolve the actual motion of individual sarcomeres within their native physiological setting— inside intact, beating cardiomyocytes—tightly integrated with quantitative modeling.

In this study, we cultured human fluorescent sarcomere-reporter hiPSC-derived cardiomyocytes (ACTN2-Citrine-CMs) on micropatterned gel substrates of different stiffnesses to investigate sarcomere dynamics under different mechanical boundary conditions. High-speed imaging, combined with AI-based analysis using our SarcAsM algorithm¹⁵, enabled us to track average and individual sarcomere motions within myofibrils with nanometer resolution.

Our analysis found strikingly complex dynamic behaviors, including rapid length oscillations, transient ‘popping’ events, and heterogeneous contraction patterns that varied both temporally and spatially along one myofibril. We show that this stochastic heterogeneity appears to arise from competitive, tug-of-war-like interactions among mechanically coupled sarcomeres within the myofibrils. To rationalize these observations, we developed a data-driven mesoscopic model which demonstrates that the emergent heterogeneity at the sarcomere-level can be explained by stochastic fluctuations in conjunction with a non-monotonic force-velocity relation with a critical force threshold. These insights into stochasticity and dynamic instability at the sarcomere level challenge traditional paradigms of myocardial contraction.

Results

Individual ACTN2-Citrine hiPSC-derived CMs on micro-patterned elastic substrates

To resolve sarcomere dynamics with high spatial and temporal resolution in PSC-derived cardiomyocytes, we used a CRISPR-engineered induced pluripotent (iPSC) stem cell line, expressing a yellow fluorescing protein as ACTN2-Citrine N-terminal fusion protein after cardiomyocyte differentiation¹⁵. Seeding of differentiated cardiomyocytes on soft elastic polyacrylamide gels functionalized with a micro-printed pattern (70 × 10 µm; Fig. 1B) promoted anisotropic myofibril assembly in uniformly elongated cardiomyocytes. Culturing cells on gels with defined elasticities (Young’s moduli: 5 kPa, 9 kPa, 15 kPa, 29 kPa, 49 kPa, 85 kPa; Table S1) allowed us to impose auxotonic loads on the cells on a scale considered to be relevant in vivo under physiological (10-20 kPa) and pathological (e.g., fibrosis; ≥30 kPa) conditions¹⁵. After a maturation period of 20-35 days, we recorded 20-30 s long movies of, in total, 1,362 spontaneously beating CMs at a frame rate of 66 Hz (Movie S1, Fig. S1).

Automated tracking of sarcomere motion at high spatial and temporal resolution

We used SarcAsM, an AI-based software tool for automated segmentation and analysis of sarcomere structure and dynamics¹⁵ for data evaluation. SarcAsM automatically identified 5,085 lines of interest (LOIs), up to 4 per cell, with well-organized registered sarcomeres (≥10 sarcomeres; Fig. 1C, S1). Along each LOI (∼800 nm in width), an intensity kymograph of Z-band motion was extracted from deep-learning (3D U-Net) processed movies, allowing for robust and time-consistent segmentation of Z-bands in noisy high-speed movies (Fig. 1C). This approach provided more accurate and robust localization and tracking of Z-band trajectories Z_i(t) of individual sarcomeres than using raw microscopy data (∼17 nm and 15 ms resolution, details see SarcAsM¹⁵). From Z_i(t), we obtained sarcomere length SL_i(t), sarcomere length change ΔSL_i(t) and sarcomere velocity V_i(t) of each sarcomere i as well as multi-sarcomere averages and for each LOI (details see SarcAsM¹⁵) (Fig. 1D-F).

Cardiomyocyte and individual sarcomere contractility as a function of substrate stiffness

The hiPSC-CMs were beating spontaneously at a frequency of 0.91 ± 0.38 Hz, which was largely independent of the substrate stiffness (Fig. S2A). Beat-to-beat intervals, though, showed increasing irregularities with increasing substrate stiffness (Fig. S2B). Contraction durations T_c shortened with increasing substrate stiffness (Fig. S2C).

In line with previous reports^16-18, total cell contraction amplitudes, quantified by the inward motion of the outermost z-bands in myofibrils, substantially decreased with increasing substrate stiffness (Fig. 1D-F). A surprising finding was that the length changes of individual sarcomeres ΔSL_i fell out of synchronization and became distinctly heterogenous, deviating strongly from the average length change with increasing substrate stiffness. Sarcomeres also showed more singular large-amplitude extensions (“sarcomere popping”) far beyond their resting lengths, often very pronounced at the end of contractions (marked with asterisks in bottom row of Fig. 1D-F and large excursions in phase-space plots in Fig. 2A-C). Despite the rapid and heterogenous motion of individual sarcomeres, the emergent contraction at the myofibril scale remained regular and smooth on all substrates.

Figure 2:

To analyze the deviations of individual sarcomere dynamics from the myofibril average during contraction cycles, we collected minima and maxima of ΔSL, i.e., maximal contraction and elongation amplitudes, and minima and maxima of V for each contraction cycle for both individual sarcomeres and averages over all sarcomeres in a given LOI (Fig. 2A-G). In total, we analyzed a data set of 5,085 LOIs recorded from 1,362 cells, each over at least 10 contraction cycles. We excluded irregularly beating (beat-to-beat variability >0.1 s) and very rapidly (>2 Hz) or slowly (<0.5 Hz) beating cells. We further limited the analysis to LOIs spanning whole myofibrils, where the average length change was strictly negative, i.e., showing contraction. Out of the full data set of 5,085 LOIs, 2,321 LOIs met these criteria. This selection ensured that we focused on rhythmically and uniformly contracting cells. As expected, maximal average contraction amplitudes , measured at the peaks of shortening in each contraction cycle, were the largest (0.16 ± 0.05 µm) on the softest (5 kPa) substrates and were close to zero (0.04 ± 0.03 µm) on the hardest substrates (85 kPa; Fig. 2D). Maximal contraction amplitudes ΔSL₋ of individual sarcomeres in each contraction cycle were also largest on 5 kPa substrates at 0.22 ± 0.09 µm and declined to 0.13 ± 0.07 µm on 85 kPa (Fig. 2E). Note that the decline of maximal individual sarcomere shortening (ΔSL₋) by ∼30% from 5 to 85 kPa was far less than the ∼75% decline of the average sarcomere contraction . In addition, while myofibrils never elongated beyond resting length, as expected for auxotonic contractions, individual sarcomeres frequently elongated well beyond their resting lengths during contraction cycles (median ΔSL₊ = 0.05 µm). The distributions of maximal sarcomere extension amplitudes (ΔSL₊) were remarkably similar between substrate conditions (Fig. 2E).

Next, we evaluated the effect of substrate elasticity on average and individual sarcomere contraction and extension velocities. The average extremal velocities and followed the same trend as the extremal average length changes and and were the largest on 5 kPa and declined steadily and strongly with increasing substrate stiffness (from 5 to 85 kPa by ∼60% for both contraction and elongation; Fig. 2F). The single sarcomere extremal contraction and extension velocities V₋ and V₊, in contrast, were far less affected by substrate elasticity and differed only by at most ∼10% (V₋) and ∼20% (V₊) for the various substrate conditions (Fig. 2G). Note that the maximal extension velocities were 29 - 38% larger than the maximal absolute contraction velocities, showing a strong asymmetry between single sarcomere contraction and extension dynamics.

High-frequency oscillatory motion of individual sarcomeres

Unlike the average sarcomere motion with clearly distinguishable contraction and relaxation phases, single sarcomeres exhibited rapid switching between slow contractile and fast extensile motions with sometimes multiple shortening and lengthening phases within one contraction cycle (Fig, 2H,I, Fig. 1D-F). To quantify this phenomenon, we analyzed the oscillation frequencies of the average and individual length changes ΔSL(t) over time, utilizing Morlet wavelet analysis, a method that makes it possible to extract the instantaneous oscillation frequencies of a signal. For average sarcomere length changes , the analysis revealed a narrow range of frequencies surrounding the beating rate in the scalogram (Fig. 2H,I). For individual sarcomeres, in contrast, we found a broader range of frequencies, with significant peaks occurring at the dominant and macroscopically visible whole-cell beating rate as well as a large secondary peak at 3 ± 0.5 Hz (Fig. 2J). This distinctive high-frequency oscillation was consistently observed across numerous LOIs and appeared to be independent of the beating rate. The intrinsic oscillatory motion of individual sarcomeres suggests a mechanism of active acto-myosin powered contraction coupled to rapid lengthening in a relaxation oscillator-like behavior¹⁹.

Correlation analysis identifies stochastic and static heterogeneity among sarcomeres

The observed heterogeneity of sarcomere dynamics, particularly on rigid substrates (Fig. 3A), might be predetermined by static non-uniformities among sarcomeres (e.g., by differences in functional myosin numbers). To distinguish this possibility from stochastic heterogeneity, we extracted the time-resolved motion ΔSL_i,k(t) of each sarcomere i during each contraction cycle k. We then quantified the similarity between these motion patterns using the Pearson correlation coefficient. For two motion patterns ΔSL_i,k(t) and ΔSL_j,l(t), this coefficient is defined as:

where t represents the discrete time index (frame number) within a contraction cycle, T is the total number of frames capturing one contraction cycle, and denotes the mean value of the sarcomere length change over the cycle.

Figure 3:

For each LOI, we calculated two correlation measures. The mutual correlation coefficient r_m measures the synchrony between sarcomeres within the same contraction cycle, calculated as r_m = ⟨r(i, j, k, k)⟩_i≠j, where the averaging is performed over all pairs of different sarcomeres (i ≠ j) within the same contraction cycle k. The serial correlation coefficient r_squantifies the consistency of individual sarcomeres across different contraction cycles, calculated as r_s= ⟨r(i, i, k, l)⟩_k≠l, where the averaging is performed over all pairs of different contraction cycles (k ≠ l) for the same sarcomere i.

The serial and mutual correlations both were maximal on 5 kPa substrates and declined steadily with increasing substrate stiffness. The decrease of mutual correlation reflects the decrease of synchrony between sarcomeres due to the tug-of-war competition imposed by the rigid constraints on the myofibril level (Fig. 3B). The decrease of serial correlation indicates an increased beat-to-beat variability of the motions of individual sarcomeres (Fig. 3B). Interestingly, the mutual correlation declined more strongly than the serial correlation, by 50% to more than 90% from 5 - 85 kPa.

We introduce the ratio R = r_m/r_s of serial to mutual correlation coefficients to distinguish static from stochastic heterogeneity (Fig. 3B,C). When R = 1, the heterogeneity is completely stochastic. Random shuffling of sarcomeres (i, j) and contraction cycles (k, l) would in that case not affect R. When R ≪ 1, the heterogeneity is largely static. In this case, mutual correlations are much smaller than serial correlations, indicating that the motion varies much more between sarcomeres than for a given sarcomere from one contraction cycle to another (Fig. 3D,E). The broad distribution of R values for all examined LOIs shows that the heterogeneity was neither fully stochastic nor fully static but rather distributed across the full spectrum from stochastic to static (Fig. 3C).

We examined R for different substrate conditions and found a strong dependence on substrate elasticity for both length change and velocity correlations (Fig. 3B). Heterogeneity was mostly stochastic for 5 kPa and 9 kPa substrates and increasingly static for stiffer substrates with non-physiological elasticities.

Sarcomere popping events are stochastically independent and not only occur in structurally “weak” sarcomeres

Beyond the heterogeneous contraction patterns, we observed discrete mechanical instabilities in the form of sudden sarcomere elongations, mostly at the end of contractions, which we termed “popping” events. Defined as sarcomere elongations exceeding 0.25 µm (Fig. 3A,F), these events likely represent passive yielding where threshold tension triggers an avalanche-like release of myosin heads while neighboring sarcomeres continue to contract or relax in a more direct manner without overshoot.

The rate of popping events showed a clear dependence on substrate stiffness (Fig. 3G). Stiffer substrates (49-85 kPa) led to higher popping rates than softer substrates (5-15 kPa), with median frequencies increasing from approximately 0.08 events per sarcomere per cycle on soft substrates to over 0.14 on the stiffest substrates. Importantly, even on the softest substrates, popping frequencies remained finite, indicating that sarcomere popping is an intrinsic feature rather than a pathological response to specific mechanical conditions. This baseline level of popping across all substrate stiffnesses demonstrates that mechanical instabilities are fundamental to sarcomere dynamics.

Popping frequency also correlated with sarcomere equilibrium length (SL₀), with longer sarcomeres exhibiting higher popping rates across all substrate conditions (Fig. 3H). However, the distribution of equilibrium lengths (gray histogram in Fig. 3H) was relatively tightly centered around 1.8 µm. Crucially, the popping rate remained consistently non-zero also across the entire range of equilibrium lengths, even in LOIs where length variability was minimal. This demonstrates that while resting length heterogeneity contributes to popping behavior, it is not the primary cause of these mechanical instabilities.

To assess the stochastic nature of popping events, we analyzed the distributions of time intervals and spatial distances between consecutive popping events. We assessed the distributions of distance (d) and time gaps (τ) for all LOIs and compared them with geometric distributions G(k) with respective event probabilities (Fig. 3I,J). Using the Kolmogorov-Smirnov test, we could not reject the hypothesis that popping events are stochastically independent for 49% of LOIs with respect to distance and 47% of LOIs with respect to time gaps (p-value = 0.01). This substantial proportion of LOIs showing stochastically independent popping events demonstrates that popping has a stochastic component, independent of deterministic factors such as sarcomere strength heterogeneity or structural variations. While the remaining LOIs show non-random clustering patterns—likely reflecting the influence of local mechanical coupling or morphological heterogeneities—the underlying stochastic mechanism appears to be intrinsic to sarcomere mechanics itself.

A mesoscopic Langevin framework for coupled sarcomere dynamics

To rationalize the observed complex sarcomere dynamics, we developed a model (Fig. 4A) describing the dynamics of each sarcomere i by an underdamped Langevin equation—a second-order stochastic differential equation. We represent a sarcomere mechanically as composed of three parallel components: an active force generator, contributing force F_a,i, a viscous damper, contributing force F_d,i, and an elastic element, contributing force F_s,i. Since we neglect physical mass, forces at all nodes between sarcomeres are instantaneously balanced, meaning every sarcomere experiences the same external load F_m. However, this external load is not instantaneously balanced by the internal forces (F_a, F_d, F_s) within each sarcomere due to the transient nature of the contraction. The resulting force imbalance distributes differently among the three elements from sarcomere to sarcomere, driving asynchronous dynamics (Fig. 4A). To capture this non-equilibrium behavior, we introduce the effective parameter μ, which acts as a quasi-inertial mass. The equation of motion reads:

where v_i and x_i are velocity and length of sarcomere is a vector of all sarcomere lengths, c(t) is the time-dependent activation, and η(t) a Gaussian noise term. Multiplying Eq. (2) by μ recovers the form of Newton’s second law (F = m · a), where represents the inertial term. While not a physical mass, μ models the characteristic time scale of the system’s kinetic response, reflecting the delayed attachment and detachment rates of actomyosin cross-bridges. This parameter defines how fast the system is drawn to the nullcline determined by the force balance, effectively smoothing instantaneous force fluctuations, representing the system’s persistence against rapid force changes.

Figure 4:

This phenomenological underdamped Langevin framework, previously successfully applied to nonlinear systems in general¹⁹ and to cell migration in particular^20,21, captures the nonlinear dynamics and emergent collective behaviors observed in our sarcomere system without explicitly modeling the underlying molecular mechanisms. This makes it particularly well-suited for studying the emergence of complex, large-scale behavior from fundamental single-sarcomere interactions.

To model the active actomyosin force F_a, we chose the 3^rd order polynomial ansatz, while the viscous drag F_d is modeled as a linear friction term:

The polynomial in Eq. (3) is a generalist functional form capable of capturing a broad spectrum of force-velocity relationships. Multiplying by the time-dependent activation c(t) makes sure that active force generation scales with calcium availability, reflecting troponin C binding and the resulting availability of myosin binding sites on actin². At high activation (c ≈ 1), F_a drives contraction, while at low activation (c ≈ 0), the viscous drag F_d dominates, governing passive relaxation dynamics..

In our model, as in previous work²², we represent the passive elastic element inside each sarcomere i using a piecewise defined function, allowing for different behaviors for shortening and lengthening. We here assume that elastic resistance during shortening is linear (provided by, e.g., microtubules²³) and quadratic during lengthening (provided by, e.g., titin²²):

Here, k₁ and k₂ are elastic constants. For simplicity and without loss of generality, we set the sarcomere equilibrium length to zero (x_eq = 0). For serially coupled sarcomeres, the total elastic force on each sarcomere is where is the external myofibril load. For our experiment of single hiPSC-derived cardiomyocytes on micropatterned elastic soft gels, we assume an auxotonic load with myofibril length change x_m = ∑_j x_j and elasticity k_l.

We model stochastic fluctuations on the sarcomere level through the random noise term η(t) = η_t(t) + η_a(t). This term combines two distinct components: thermal noise η_t(t) and active noise η_a(t), originating from the stochastic dynamics of the collective action of the myosin motors in the sarcomere. As shown previously²⁴, these fluctuations arising from the stochastic motor kinetics can be approximated as delta-correlated:

with system-specific magnitude η_a,0, the fraction of active motors n and the duty ratio of motors r. These fluctuations originate from the discrete, stochastic switching of individual motors between bound and unbound states²⁴. For the active fluctuations, we model n with a sigmoidal function with steepness β dependent on velocity v dependent on activation c(t), assuming active shortening and passive elongation of sarcomeres:

Data-driven parameter optimization reproduces complex sarcomere dynamics

Our model comprises 10 parameters. To specify the parameters, we simulated 20 serially connected sarcomeres over 10 contraction cycles using the Euler-Maruyama scheme²⁵. The cyclic activation was modeled using a simplified custom activation function c(t) based on a clipped sine with a beating rate matching the experimental data (Fig. 4B). Since analytical fitting of such non-linear, stochastic and time-dependent systems to experimental data is impossible, we employed Differential Evolution²⁶, a robust gradient-free optimization method, to identify the physical parameters generating the experimentally observed dynamics. The objective was to minimize a loss function based on the Kolmogorov-Smirnov (KS) coefficient, which quantifies the similarity between the distributions of simulated and experimental data on a scale from 0 (identical) to 1 (no overlap).

The model evaluation compared the distributions of sarcomere length changes and velocities from simulations with representative experimental LOIs from substrates (5, 15, and 85 kPa) covering the full range of mechanical loads. To account for the transient activation dependency, we divided the trajectories into four distinct phases: quiescent, early activation, mid-activation, and late activation, and compared the distributions for each phase using the KS metric. To match the experimental temporal resolution, model data was subsampled to 16 ms frame intervals. For each parameter set during optimization, all three substrate conditions were simulated, and the KS coefficient was computed for each condition and phase. The final loss was then calculated as the mean loss over all three conditions and four contraction phases, yielding a single metric quantifying the overall model-experiment agreement. The 10 model parameters were optimized within prespecified bounds (Table S2).

This optimization process yielded a model in excellent agreement with the experimental data (Fig. 4D-F). The average KS coefficient across all conditions and phases was ∼0.05, confirming a high degree of similarity between the simulated and measured distributions of individual sarcomere lengths. Ultimately, we inferred a single, unified model capable of replicating dynamics across all substrate elasticities by adjusting only the substrate stiffness parameter. The model correctly captures the extremal trends (minima and maxima) of decreasing sarcomere displacement and velocity with increasing substrate stiffness at both the individual and average sarcomere levels (Fig. 4G, H). While the model predicts slightly stronger average contraction and velocity extrema than measured, the stiffness dependence of these extrema is consistent with experimental observations (gray data in Fig. 4G, H). Note that the model correctly replicates the emergence of smooth and periodic contraction at the myofibril level from the highly stochastic and heterogeneous dynamics of individual sarcomeres (Fig. 4D-F).

Non-monotonic force-velocity dynamics drive tug-of-war and limit cycles in coupled sarcomeres

Our data-driven inference procedure reveals, without making this an a priori assumption, that the active force F_a exhibits a non-monotonic force-velocity relationship (Fig. 4C, Movie S2), with branches separated by an unstable region of negative slope (dF/dv < 0), corresponding to negative friction. This relationship is dynamically modulated by the activation c(t), creating a complex time-dependent force landscape that governs sarcomere behavior.

The model is effectively analogous to a Van-der-Pol relaxation oscillator in force-velocity space¹⁹. As total myofibril shortening proceeds, the external mechanical load from the substrate increases, slowing down contraction velocities of all sarcomeres. The non-monotonic force-velocity relationship implies that there is an intrinsic threshold force at which sarcomeres become dynamically unstable. Physically, this corresponds to a tipping point where the stochastic unbinding of a few myosin cross-bridges triggers an avalanche-like release of the remaining more strongly loaded cross-bridges, switching the sarcomere from contraction to rapid elongation.

This instability has different consequences in the different phases of the contraction cycle: During the initial contraction phase, joint sarcomere shortening drives the velocity toward the instability threshold. Through dynamic selection, stochastic fluctuations (modeled by η) then push random sarcomeres past the tipping point. Because activation remains high, this elongation is transient; as a particular sarcomere lengthens, load quickly redistributes to its passive elements and while other sarcomeres take up the slack, it ends up back on the contracting branch. This rapid switching generates the stable limit-cycle oscillations observed in individual sarcomeres (Fig. 2H-J).

During relaxation (declining activation), the instability threshold decreases (Fig. 4C). On stiff substrates, where external load remains high, this lower threshold causes a larger fraction of sarcomeres to slip past the instability point. Unlike in the first regime, the declining active force allows for large, sustained “popping” excursions (Fig. 4D-F). We observed that only a subset of sarcomeres pop (∼10%), which effectively buffers the load for the remaining sarcomeres, allowing them to relax smoothly without over-extension.

Stochastic fluctuations η(t) (Eq. 6) are essential to reproduce the heterogeneous sarcomere behavior seen experimentally, including beat-to-beat variability, competition, and popping events. Removing the stochastic drive from the equation of motion (Eq. 2) in simulations collapses trajectories, inconsistent with experiments (Fig. 5A,B). To probe robustness against fixed force non-uniformities, we introduced static variability by scaled each sarcomere’s active force F_a with a constant, factor drawn for each sarcomere in the chain of 20 from a distribution with standard deviation σ_a. Functionally, this assigns each sarcomere a slightly different maximum force threshold. This force variability alone biases trajectories and creates deterministic static heterogeneity during contraction according to each sarcomere’s strength (Fig. 5C), also inconsistent with our data. Adding stochastic fluctuations leads to trajectories most closely resembling the experimental data (Fig. 5D).

Figure 5:

Correlation analysis (Fig 5E; compare with experimental data in Fig. 3B) of the simulated data with different levels of static heterogeneity σ_a shows, as expected, that without stochastic fluctuations, the serial correlation of sarcomeres is always 1, meaning that individual sarcomeres always behave the same way through contraction cycles. In contrast, simulations with stochastic fluctuations maintain stochastic heterogeneity (R > 0.8) even with substantial static force heterogeneity (up to σ_a ≈ 0.2) (Fig. 5E).

Discussion

We combined high-speed imaging of human induced pluripotent stem cell (iPSC)-derived cardiomyocytes with detailed AI-based quantitative analysis. We discovered rich subcellular dynamics that remain hidden at the whole-cell level. While overall cell contraction cycles were regular, individual sarcomeres displayed rapid oscillations, reversible popping events, and stochastic beat-to-beat variability within the same myofibril. Increased substrate stiffness reduced overall cell contraction amplitudes, but had only a weak effect on single-sarcomere length-changes and contraction velocities. Rigid substrates instead promoted tug-of-war-like interactions and loss of synchrony between sarcomeres along myofibrils.

To interpret these observations mechanistically, we developed a data-driven, mesoscopic modeling framework describing transient, non-equilibrium dynamics and mechanical coupling among serially connected sarcomeres. The model relies on two components that accurately reproduce the experimentally observed dynamics: a non-monotonic force-velocity relationship that generates a dynamic instability at a critical force and intrinsic stochastic force fluctuations at the sarcomere level. This approach differs from classical modeling of skeletal muscle based on a force-length relationship. We find that our model is more appropriate than such models for cardiomyocytes that operate in a regime of rapid transient activation-relaxation cycles, viscous drag, and velocity-dependent force generation.

The dynamic instability central to our framework arises from the non-monotonic shape of the force-velocity relationship. Note that this functional form was not an initial assumption but that it emerged when optimizing the parameters of the very general modeling ansatz to fit our experimental data. This result connects to foundational, reductionist theories of coupled molecular motors working against an external load^27-30, which predict that at a critical load, a collective, avalanche-like unbinding of motors from their track (e.g. myosin heads from actin in a sarcomere) induces a dynamic phase transition from active forward motion against the load to rapid passive slipping. Such instabilities have so far only been reported in reconstituted motility assays^29,31. We have here observed an analogous instability in cardiomyocytes manifesting itself in the observed switching between contraction and expansion and the popping events of individual sarcomeres. If sarcomeres followed a monotonic force-velocity relationship, increasing load would simply cause uniform slowing or stalling. Small non-uniformities caused by fluctuations would decay, returning the system to a uniform state. The non-monotonic force-velocity relationship, instead, amplifies such fluctuations, driving the system into limit cycles that produce the high-frequency oscillations and stochastic heterogeneity we observe (Fig. 2G).

What could be the physiological roles of the observed phenomena? The popping events are the most dramatic outcome of the instability-amplified stochastic noise. Their reversible nature, ubiquity across all tested conditions, and prevalence at the end of contractions - where declining activation reduces the fraction of attached cross-bridges - are all consistent with a dynamic instability rather than structural failure. We therefore view popping as an intrinsic outcome of transient, non-equilibrium sarcomere dynamics. While possibly similar popping has been described in skeletal muscle³², particularly in the context of residual force enhancement³³, the observation in rhythmically contracting cardiomyocytes is a key new finding of this study. A physiological function of popping might be that intermittent yielding in a subset of sarcomeres accelerates myofibril-level relaxation toward the onset of the next diastole, helping tension decay even while individual sarcomere lengths remain non-uniform. A physiological role of stochastic heterogeneity of yielding might be the limitation of sarcomere damage. Our data indicate a mixture of static and stochastic heterogeneity. Some sarcomeres exhibit persistent differences suggestive of structural or compositional non-uniformities, yet many myofibrils show strong, cycle-to-cycle stochastic variability, including popping. As our simulations demonstrate, even in the presence of low to moderate sarcomere strength variability, the interplay between stochastic fluctuations and the dynamic instability prevents deterministic yielding. The interplay between intrinsic fluctuations and the dynamic instability randomizes yielding events across cycles and sarcomeres, preventing the repetitive overload of any single element and reducing the likelihood of localized damage.

Stochastic heterogeneity thus provides a robust mechanism for maintaining structural and functional homeostasis. Healthy adult cardiomyocytes—with their more ordered myofibril and sarcomere architecture compared with stem-cell-derived cardiomyocytes³⁴—are expected to show smaller static non-uniformities while still exhibiting stochastic heterogeneity. While speculative for human cardiomyocytes, this hypothesis is supported by in-situ measurements of single-sarcomere dynamics in mouse myocardium, which found stochastic heterogeneity⁶.

In diseased and disordered myofibrils, static non-uniformities may overwhelm the protective stochastic mechanism ^14,34,35 and instability may then be consistently triggered at the most vulnerable sarcomeres, concentrating mechanical stress and leading to fatigue and maladaptive remodeling. This perspective provides a more nuanced view than those based on fixed “weak-strong” contrasts¹⁴, suggesting that pathology may arise from local structural weakness, but that that weakness would have to overwhelm the protective load-sharing mechanism. The phenomenon that amplified random noise can enhance responsiveness and stabilize dynamics has also been described in systems showing stochastic resonance³⁶.

The presented framework is mesoscopic by design, omitting explicit molecular details such as calcium kinetics or individual cross-bridge states. Instead, we adopted a top-down approach, inferring the emergent mesoscopic physics directly from the experimental data. We utilized a mesoscopic Langevin framework in the hydrodynamic limit, where the deterministic force term represents the collective average of the underlying molecular kinetics and mechanical interactions, while the stochastic term accounts for unresolved microscopic fluctuations, following the general principles of active matter theory applied to actomyosin systems^37,38. The strength of this coarse-grained approach is that it reveals emergent physical principles—such as the non-monotonic force-velocity relationship—that would remain obscured in parameter-heavy molecular models. The model’s predictive power, with only 10 parameters capturing dynamics across three substrate stiffnesses and four contraction phases, validates this strategy.

Our combined experimental and modeling platform offers a powerful tool to dissect the specific molecular drivers of collective sarcomere dynamics in future studies. By applying pharmacological perturbations or specific mutations, one could link distinct molecular defects — such as altered cross-bridge kinetics, calcium sensitivity, or titin stiffness — to quantitative changes in the model’s physical parameters. This would establish a rigorous mapping between molecular function and myofibril stability, enabling a systematic analysis of how specific pathologies disrupt the contractile homeostasis of the heart.

Materials and Methods

Generation and culturing of hiPSC ACTN2-Citrine-derived cardiomyocytes

Cardiomyocyte differentiation of a ACTN2-Citrine reporter hiPSC line³⁹ was performed according to Tiburcy et al.⁴⁰. hiPSC-ACTN2-Citrine-derived cardiomyocytes were cultured in 6 well plates (Cat 3516, Corning) in serum-free “cardio” medium (0.4 mM Ca²⁺, RPMI 1640 with GlutaMAX (Cat 61870, Invitrogen), 1% Penicillin/Streptomycin (Cat 15140, Invitrogen), 2% B27 supplement (Cat 17504-044, Invitrogen) at 37 ºC in a 5% CO₂ incubator with culture medium changes at every other day. For re-seeding, cells were detached using Accutase^® digestion medium (StemPro^® Accutase^® cell dissociation reagent (Cat A11105-01, Gibco), 0.025% Trypsin (Cat 15090-046, Gibco), 20 µg/mL DNaseI (Cat 260913, Calbiochem) for 15-20 min at 37ºC. Digestion was stopped using “cardio” medium supplemented with 5 µmol/L Rock Inhibitor (Stemolecule Y27632, Cat 04-0012-10, Reprocell) three times the volume of the Accutase^® mix. Cell clumps were separated using a 100 µm cell strainer. Cells were seeded using ∼150,000 cells per micropatterned patterned substrate of ∼1 cm² size and cultured 24 hours in “cardio” medium with 5 µmol/L Rock Inhibitor. Cardiomyocytes on soft gels were maintained in serum-free “cardio” medium with daily medium changes at 37 ºC in a 5% CO₂ incubator for up to 30 days. Cardiomyocytes were imaged after a maturation period of 20-30 days post seeding on the soft gels.

Sample preparation and measurement

For live-cell video imaging, the soft substrates were mounted in a custom-built holder for round ∅25 mm glass cover slides with serum-free medium. Movies of beating cardiomyocytes were obtained with a confocal microscope (TCS SP5 II, Leica, Germany) at 37ºC and 5% CO₂. We used an 8,000 Hz resonant scanner with bidirectional scanning mode and recorded up to 20-30 s long movies of 1024 × 200 pixels at 67 frames per second, i.e. a temporal resolution of 15 ms.

Fabrication of cell-adhesive micropatterns on polyacrylamide soft gels

Silicon wafers (Microchemicals GmbH, Ulm, Germany) were coated with SU-8 photoresist (Series 3005, MicroChem, Newton, USA) using a two-step spin-coating process. The wafers were then exposed to UV light through custom photomasks (Compugraphics, Jena, Germany) and developed to create photoresist masters (details in Ref³⁹). Polydimethylsiloxane (PDMS) stamps were produced by mixing PDMS and curing agent 10:1 (Sylgard 184 kit, Dow Corning) and pouring it onto the photoresist masters. After degassing and curing, the PDMS was cut and peeled off to create the stamps. Micropatterned polyacrylamide gels were prepared by treating PDMS stamps with plasma to make them hydrophilic, then incubating them with 0.1 mg/ml Synthemax™ (Cat 3535, Corning). The protein-coated stamps were placed on plasma-cleaned glass coverslips and weighted to transfer the protein. Gel solutions of acrylamide and bis-acrylamide were prepared to achieve different elastic moduli (Table S1), measured with a rheometer (Physica MCR 501 Rheometer, Anton Paar, Austria). The gel solution was polymerized with the Synthemax-patterned glass on top. The gels were stored in PBS, then washed before use.

Microscopy data analysis and statistical analysis

For processing of the live-cell confocal movies, and tracking and analysis of sarcomere trajectories, our Python package SarcAsM (Sarcomere Analysis Multitool) was used¹⁵. All additional analyses and illustrations were created using custom scripts written in Python. All data is displayed as mean ± standard deviation unless indicated otherwise. Whenever applicable, we employed the Kruskal-Wallis test, a non-parametric method, for hypothesis testing, and proceeded to a post-hoc multi-comparison using Dunn’s test, setting the threshold for statistical significance at a p-value of less than 0.01. All condition differences were deemed significant except where explicitly marked as not significant (n.s.).

Computational model

We implemented the mesoscopic model of coupled sarcomeres in Python 3.12, making use of the NumPy⁴¹ and SciPy⁴² libraries. The model describes sarcomere dynamics through serially coupled underdamped Langevin equations, which is an approach suited to capturing emergent collective behaviors at the myofibril scale without explicitly modelling molecular-level details. We solved the system of stochastic differential equations numerically using the Euler-Maruyama scheme²⁵ with a time step of 2 ms. The dimensionless load parameter k_m was provided as a simulation input, defined for simplicity as 0.1 times the substrate elasticity (in kPa). The remaining model parameters were estimated using Differential Evolution²⁶, a global optimization algorithm, to minimize the Kolmogorov-Smirnov distance between simulated and experimental distributions.

Data availability

All data, analysis scripts, and source code supporting this study are publicly available. The primary software packages developed for this research, SarcAsM and the computational model, are available on GitHub at https://github.com/danihae/SarcAsM and https://github.com/danihae/SarcomereModel, respectively. The complete experimental dataset, alongside code to reproduce the analyses in this manuscript, is archived on Zenodo (https://doi.org/10.5281/zenodo.17564384). This repository includes a representative dataset of 25 raw microscopy movies with their corresponding processed line-of-interest (LOI) sarcomere trajectories for each of the six substrate stiffness conditions, as well as the aggregated data frame used to generate all figures.

Acknowledgements

We thank Florian Rehfeldt, Andrei Vilfan, Lev Truskinowsky, David Brückner, Chase Broedersz and Pierre Ronceray for helpful discussions. We gratefully acknowledge the use of the microscopy facility of the Max Planck Institute for Multidisciplinary Sciences for access to cell culture and high-speed confocal microscopy. CFS and DH would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program “New statistical physics in living matter: non equilibrium states under adaptive control”. DH acknowledges the support from the German Academic Foundation (Studienstiftung des Deutschen Volkes) for providing a doctoral fellowship and the Campus Institute for Data Science (CIDAS) at the University of Göttingen for awarding a postdoctoral fellowship. The research of CFS was supported by the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) with the ERC grant agreement nº340528. WHZ acknowledges the support from the DZHK (German Center for Cardiovascular Research), the German Federal Ministry of Education and Research (IndiHEART; 161L0250A), the German Research Foundation (DFG SFB 1002 C04/S01, IRTG 1816, RTG 2824, EXC 2067-1), and the Fondation Leducq (20CVD04).

Additional information

Author contributions

D.H., T.D., W.H.Z. and C.F.S. conceived and designed the study. D.H., L.H., and T.D. developed the methodology. W.H.Z. and C.F.S. supplied the essential resources. D.H., L.H., and T.D. performed micropatterning and cardiomyocyte experiments. D.H. and K.N. analyzed experimental data. D.H. conceived and implemented the theoretical model. D.H. created the visualizations. The original draft was written by D.H., and all authors contributed to the review and editing of the manuscript.

Funding

EC | European Research Council (ERC) (340528)

• Christoph F Schmidt

German Center for Cardiovascular Research

• Wolfram Zimmermann

German Academic Foundation

• Daniel Haertter

Campus Institute for Data Science

• Daniel Haertter

German Federal Ministry of Education and Research (IndiHEART; 161L0250A)

• Wolfram Zimmermann

German Research Foundation (DFG SFB 1002 C04/S01)

• Wolfram Zimmermann

German Research Foundation (IRTG 1816)

• Wolfram Zimmermann

German Research Foundation (RTG 2824)

• Wolfram Zimmermann

German Research Foundation (EXC 2067-1)

• Wolfram Zimmermann

Fondation Leducq (Leducq Foundation) (20CVD04)

• Wolfram Zimmermann

Additional files

Supplemental Information

Movie1

Movie2