Theory of frequency response of mechanically driven cardiomyocytes

We theoretically predict and compare with experiments, transitions from spontaneous beating to dynamical entrainment of cardiomyocytes induced by an oscillating, external mechanical probe. In accord with recent experiments, we predict the dynamical behavior as a function of the probe amplitude and frequency. The theory is based on a phenomenological model for a non-linear oscillator, motivated by acto-myosin contractility. The generic behavior is independent of the detailed, molecular origins of the dynamics and, consistent with experiment, we find three regimes: spontaneous beating with the natural frequency of the cell, entrained beating with the frequency of the probe, and a “bursting” regime where the two frequencies alternate in time. We quantitatively predict the properties of the “bursting” regime as a function of the amplitude and frequency of the probe. Furthermore, we examine the pacing process in the presence of weak noise and explain how this might relate to cardiomyocyte physiology.


A Biophysical origin and dynamics of spontaneous beating
We start with the generic equation for the non-linear Rayleigh oscillator [1]: We discuss here three possible scenarios that motivate the biological origin of this equation for the oscillating degree of freedom χ(t): 1. χ(t) represents the physical displacement of the contractile machinery (sarcomeres) coupled to the deformations induced by the probe (f ext (t)). This model was discussed in detail by Jülicher and Prost [2,3] in the context of acto-myosin displacement within the sarcomere. The elastic response of the acto-myosin bundle is manifested in the spring constant k. The first order friction coefficient, ρ accounts for both dissipation of energy via the viscous environment as well as for the input of energy by the ATP dependent activity of myosin. Therefore, unlike purely damped oscillators, ρ can have a negative sign (when activity dominates friction) which allows for spontaneous contraction oscillations with frequency ω c = k/γ. If only linear terms are kept, energy input (the case of r < 0) would result in a diverging amplitude. To maintain stability, the overall energy input must dissipated by a higher order friction term λ > 0, and is related by Jülicher and Prost to the loading rate of myosin heads. Lastly, f ext (t) is the mechanical force applied to the contractile unit via the cellular adhesions that connect the sarcomeres to the substrate, such as focal adhesions at the cell edges [4] or integrin-mediated adhesions (by costameres) along the line of sarcomeres [5].
2. χ(t) represents the ionic calcium concentration in the cytosol, and the external force applies tension to the either the sarcolemma, or directly to the sarcoplasmic reticulum (SR) -affecting calcium release. In this context, ρ relates to auto-catalytic activation of ryanodine receptors (RyR) on the SR, and the higher order restoring term λ > 0 represents deactivation of RyR (which is still not fully understood) once a certain threshold of cytoplasmic concentration is reached [6,7]. The spring constant k entails numerous mechanisms such as SR Ca 2 +-ATPase activity, sarcolemma Na + -Ca 2 + pumps and mitochondrial Ca 2 + uniports, all of which working to reduce calcium concentration to baseline levels. Lastly, the external force f ext is coupled to the sarcolemma (or the SR), through integrin adhesions [5]. The tension applied to the sarcolemma (or the SR membrane) modulates the activity of mechano-sensitive protein complexes that modulate calcium release from SR [8].
3. χ(t) represents the calcium concentration as in the second scenario, but the external force is coupled via the adhesions of the substrate to sarcomeres. In this case, the deformations in the substrate by the probe modulate the binding kinetics of actin-myosin, which, in turn, either apply further tension to the stress sensitive sarcolemma or SR, or reduce the amount of calcium bound to troponin at rest. This modulates the effective ionic concentration, and further paces the calcium concentration χ(t). In this scenario, contractility is a crucial intermediate step in entrainment of calcium.
For each of these cases, the acceleration term Eq. A.1, proportional to γ, is not at all related to inertia (mass) because both the cell interior and exterior are aqueous solutions where frictional dissipation is dominant and inertia is negligible [9,10]. Rather, this term originates in the dynamical active process of binding and unbinding of myosin motors to actin filaments (first scenario), or from the dynamics of opening and closing of RyR channels (second and third scenarios). The separation of the time scales associated with binding/unbinding of myosin, or opening/closing event of RyR, and the time scale of apparent oscillations yields an "effective mass" γ > 0 term that is a function of the different kinetic rate constants in the system (as shown in detail in Ref. [3]).
To demonstrate this, we study a simple system that couples contraction of a sarcomere (x(t)) to the actin bound myosin concentration (µ(t)) (this is analogous to the derivation of Jülicher and Prost [2,3], but perhaps more transparent):μ Here, κ + /κ − are binding unbinding rates of myosin to actin, µ tot is the overall myosin concentration available for binding in a sarcomere, η f is the friction coefficient and k is the elastic response of the sarcomere modeled as a damped spring. A non-linear friction term proportional to Λ > 0 is introduced, for reasons explained below. In this simple picture, the contraction increases the overlap region between myosin and actin, allowing for more myosin heads to bind feeding back into the myosin rate equation (as expressed by the term proportional to a 1 ). The force generated by the sarcomere is proportional to the number of bound myosin heads (hence a 2 ). The solution for Eq. A.2 for the myosin concentration, as a function of the general, time dependent displacement, is given in terms of the Green's function G(t): is the solution of the homogeneous equation and µ 0 some arbitrary initial condition.
The myosin binding and unbinding rates are typically much faster (∼ 0.01 − 0.1 Hz) [11] than the time-scale of associated with oscillations (ω c ∼ 1 Hz). In the context of the above model, this implies that there are numerous binding/unbinding events of myosin heads in one contraction cycle of the sarcomeric unit. The Green's function for Eq. A.4 is given by G(t−t ) = exp(−(κ + +κ − )(t− t )) which is sharply peaked around t = t (since k + , k − << ω c , the time scale for changes in displacement). If the Green's function is sharply peaked about the time t = t, we can expand the slowly varying velocityẋ(t ) about t < t, and introduce it back into Eq. A.4: The feedback between the velocityẋ(t) and the position x(t) is proportional to the integral of the Green's function (zeroth order moment). The inertial term is associated with the first moment of the Green's function, which accounts for force accumulation due to the positive feedback between contraction and the increase of bound myosin. In our simple picture, the Green's function is , which is indeed sharply peaked around t = t. Using this result, the equation for myosin concentration in our example scales as: which, when introduced back into Eq. A.3 for the displacement yields: With γ * = a 1 /(κ + + κ − ) 2 and r * = (η f − a 1 /(κ + + κ − )). Note that the activity of myosin results in an inertial term with an effective mass γ * . The "mass" γ * is related to the cooperative dynamics of myosin binding and the active force they apply to contract the sarcomere. Note that r * can be negative if the input of energy to the oscillator from the activity dominates over friction (i.e. -η f < a 1 /(κ + + κ − )). In that case, the sarcomere will spontaneously contract. Note that the non-linear frictional dissipation term (Λ > 0) must be present since otherwise the input of energy for the case where r * < 0 will causes diverging oscillation amplitudes.
With the inertial term derived, and a higher order dissipation term Λ > 0 introduced to suppress unbounded growth, the system displays spontaneous oscillations for values r < 0. This is inherently different from a regular damped oscillator (r > 0), where any dissipation causes any initial displacement or velocity to eventually decay to 0. In the absence of external pacing f ext (t) = 0, for the case where the spring constant is large (stiff spring) relative to the dissipative terms (r, Λ k), one solves Eq. A.1 and finds that the system oscillates with the natural frequency ω c = k/γ. Therefore, we insert the solution χ(t) = A(t) cos(ω c t) into Eq. A.1 to find: For long times, we expect the amplitude to relax to its steady-state value, A s . The steady-state is defined by the oscillator state when the amplitude does not change in time. We therefore setȦ andÄ equal to zero, and collect terms of sine and cosine: For Eq.A.9 to be valid for all times, the terms multiplying the sine and cosine separately must equate to zero: The first condition determines the oscillation frequency in steady-state, and the second determines the average steady state amplitude A s = 4r/3Λω 2 c as shown in the paper.

B Dynamics of a paced cell B.1 Derivation of the Adler equation
We now introduce periodic forcing by an external probe/cell f ext (t) = F cos(ω p t), which oscillates at a frequency that is, in general, different from the spontaneous beating frequency ω c . Rescaling all the coefficients of Eq. A.1 by the effective mass γ, so that ρ = r/γ, λ = Λ/γ, A p = F/γ (same notation as in the main text) we get:χ For the entrained regime, where the system oscillates at the probe frequency ω p , we write the solution in complex form: are respectively, the amplitude and phase of the oscillating displacement χ(t).
Since the observed time that it takes a spontaneously beating cell to synchronize with the probe is much longer than the characteristic beating time (15 minutes compared to a time scale of 1 second, corresponding to ∼ 1 Hz oscillations) we use the method of averaging -where we average over the slowest mode of oscillation [12]. For an arbitrary function of time g(t) with period T , this is written as: Inserting Eq. B.12 into Eq. B.11, and averaging over the slowest mode (with a period 2π/ω p ) we write the dynamics of a(t) as: We now define the time scale a/ȧ ∼ τ which is much longer than the time scale for oscillations 1/ω p (since the amplitude and phase are slowly varying). Because this time scale is long, we consider only terms up to order 1/τ in Eq.B.14. All terms higher order thanȧ/a are neglected in the long-time limit. Using the definitions in Eq. B.12 and comparing real and imaginary terms, we obtain two dynamical equations for the phase amplitude A(t) and phase φ p (t): For weak forcing A p , the dynamical equation for the amplitude becomes independent of phase, and the steady-state amplitude has a similar form to that of a spontaneously beating cell, i.e. -A s = 4ρ/3λω 2 p (where ω c is replaced by ω p due to the change in frame of reference).
We scale the probe amplitude (A p ) by the amplitude of spontaneous oscillations (A s ) and define the scaled probe amplitude f s = A p /(2ω p A s ) as in the main text. Focusing on the case where the spontaneous frequency and the probe frequency are not very different (i.e. -ω c ∼ ω p ) we obtain the Adler equation (Eq. 2) as presented in the text. Note that in the experiments performed by Nitsan et al. [13] the amplitude ratio (equivalent to f s ) was kept fixed, so that the beating amplitude of the cell and the beating amplitude generated by the probe in the vicinity of the cell were roughly equal (A s ∼ A p ). According to our theory, this implies that the transition from entrained beating to bursting should occur when the frequency of the probe is roughly twice the frequency of the cell, the same order of magnitude as observed by Nitsan et al. [13].

B.2 Solving the Adler equation
First, we rescale time by the detuning ∆ω = (ω c − ω p ), so that t * ≡ ∆ωt. This allows us to write the dimensionless form of Eq. 2: where we define the coupling strength Q = f s /∆ω as in the main text. This dynamical equation has two bifurcation points, Q cr = ±1 where the steady-state solution change their character. For |Q| ≥ 1 the equation has a steady-state solution corresponding to φ ss = arccos(Q −1 ). This describes the regime where cell beating is entrained to the probe frequency ω p . However, for |Q| ≤ 1 there is no steady state solution, and the phase grows with time. This corresponds to a beating cell not fully entrained to the probe. In the limit of very weak forcing, or a large frequency difference ∆ω, the phase grows linearly in time, which means that the cell beats with its own spontaneous frequency ω c . However, for a value of Q close to the bifurcation point (|Q| → 1) we expect the cell to be entrained most of the time, with occasional periods of non-entrained behavior (which, for now, we assume nothing about).
To find an expression for the phase φ p , we integrate Eq. B.16 for the two regimes of Q. This yields: where we define the characteristic time τ = 2π/(∆ω (1 − Q 2 )).For |Q| < 1, the square roots are imaginary which transforms the hyperbolic tangent to a regular one. Therefore, the equation becomes periodic in time with m = 1, 2, 3... the period index, defined for the time intervals πm − τ /2 < t < πm + τ /2. Note that the time scale τ cannot explain the observed time of about 15 minutes required to entrain spontaneously beating cells, as this would imply Q extremely close to the critical value Q cr = ±1, which is not necessarily the case in the experiments. In the main text, we discuss a more robust scenario for the long time needed to transition from spontaneous to entrained beating.
We now derive the time regime for Q < 1 for which the cell is approximately entrained and the time regime for which the cell beats spontaneously (with ω c ). To examine the periodic behavior of the solution in Eq. B.17 we arbitrarily set φ 0 = 0 without loss of generality, and evaluate the time derivative of Eq. B.17.
For Q → 0 Eq. B.18 goes to unity, which is the limit where the cells beat with their spontaneous frequency. As Q is gradually increased, the phase shows distinct time regimes of intermittent slow and rapid changes (see Fig. 1). We show in the main text that the time regimes where the derivativeφ p (t) is nearly zero corresponds to an entrained cell (beating with ω p ) while the time regimes whereφ p (t) ≥ 1 corresponds to a non-entrained cell (beating with ω c ). To estimate the fraction of each period for which the cell is approximately entrained (the temporal regime whereφ p ≈ 0), we expand Eq. B.18 in time t * , around any of the mid-points (see arrow in Fig. 1C) of the "slow" change in phase: Thus in the vicinity of the mid-point, the change in phase is parabolic, and an estimate of the point of crossover (time t * c ) fromφ ≈ 0 toφ ≈ 1 occurs when the second term in Eq. B.19 approximately equals the first. This yields t * c = 2/ (Q(1 − Q)), which for Q → 1 goes to infinity; namely, the cell becomes entrained for the entire time of its beating cycle. In the opposite limit (small Q), we expand around the region whereφ p ≥ 1 of any one cycle to find: Using an analysis similar to the one for Q → 1, we see that when Q → 0 almost all of the beating cycle is characterized byφ ≈ 1 (beating with spontaneous frequency). By comparing the first and second term as before, we get t * c = 2/(Q(1 + Q)) which diverges for Q → 0. This means that in the limit of very weak pacing force, the fraction of the beating cycle for whichφ p ≈ 1 is almost unity and the duration for which the cell is entrained is (φ p ≈ 0).

C Dynamics of adaptive response
The experiments show that the time required for a cell to transition from spontaneous to entrained beating is about 15 minutes [13], which is two to three orders of magnitude longer than the duration of a single beating cycle (∼1 sec). The observed long time for the transition from spontaneous beating may be due to biological adaptation needed to transmit mechanical information (probe-induced strain of the cellular adhesions) to the cytoskeleton, resulting in a change in the beating frequency. We therefore apply a phenomenological description since the biological details are not known and are beyond the scope of the theory presented in this paper. Below, we consider the case where the adaptation process is manifested in the ability of the cell to sense the periodic deformations by the nearby probe.
We start by modifying the pacing force in Eq. B.11 to include exponential dependence in time, with a characteristic time scale τ a (related to the biological modifications within the cell): where τ a is the adaptation time which experiments suggest, is much longer than that of a single beating cycle (i.e.∼ 1/ω c , 1/ω p ). Anticipating that the steady-state will be entrained beating, we use a reference state where the phase φ p (t) is relative to ω p t and write χ(t) = A(t) cos(ω p t + φ p (t)), with A(t) and φ p (t) the time dependent amplitude and phase respectively. Following the same derivation as before (phase slowly varying, stiff spring compared to dissipative response) we average over a single cycle and separate real and imaginary parts as in Eq. B.14. This yields two dynamical equations for the phase and amplitude: where the (time dependent) correction, due to the finite adaptation time τ a are: Note that if we take ω p τ a 1 (since experiments suggest that the adaptation time is of the order of minutes, and the frequency is of the order of seconds), we can expand the corrections in the small parameter (ω p τ a ) −1 and simplify: For small values of probe amplitude A p , the steady-state amplitude of the paced cell does not vary much from that of the spontaneously beating cell. If the amplitude reaches its steady-state value faster than the phase, the two equations of Eq. C.22 are decoupled, and we are left with a single dynamical equation for the phase: From Eq. C.24 we see that when the adaptation time is much longer than the beating time ω p τ a 1, the function ψ(t) drops to zero (from its initial value of π/2 at t = 0 in a time of order ω p which is much smaller than the adaptation time. We thus consider times greater than this value and set φ p (t) to be zero in Eq. C.25. For simplicity of notation, we consider the case of small frequency differences and write: which is the modified Adler equation presented in the paper. One can immediately see that for short times t τ a the pacing is effectively zero, and for long times t τ a the pacing is maximal.
In Fig 2 we plot the time evolution of the phase φ p and its time derivativė φ p (both normalized for convenience) as a function of reduced time t * = ∆ωt, for various values of the scaled adaptation time ∆ωτ a . One can see that as τ a increase, the phase relaxes to a constant value at longer and longer times. This means that the time required to achieve synchronization becomes much longer. Also, note that in all cases where the cell is eventually synchronized, there is a slight overshoot (corresponding to a single root for the matching derivative) followed by a smooth relaxation to the steady-state synchronized phase. This peak time marks the transition from oscillatory, to relaxation behavior. We therefore use the peak time t c as a measure of the time required to synchronize the cell.
In Eq. 2, the transition from oscillatory to relaxation dynamics happens when the parameter Q changes sign. This suggests that for the analogous Eq. C.26, the transition would occur at time t tr for which: Eq. C.27 can be solved analytically for all Q > 1, which yields: where W (z) is the Lambert (Product-Log) function. In Fig. 3A we show the transition time t c estimated numerically, plotted against the (corrected) transition time t tr evaluated analytically from Eq. C.28, for varying values of Q and τ a . As predicted by our analysis, numerical solutions of the equation collapse rather well onto a linear curve of slope 1. In Fig. 3B we plot the ratio of the correction w(Q, τ a ) and the scaling argument given in the main text (derived neglecting w (Q, τ a )). For Q ∼ 1 the correction is negligible (the ratio is 0), showing that the scaling argument is valid. As Q increase, w(Q, τ a ) becomes more and more negative, effectively shortening the transition time. For Q 1 the scaling argument and w(Q, τ a ) are comparable in size and opposite in sign, leading to a shortened transition time that scales as ∼ τ a /Q.

D Dynamics of stochastic beating D.1 Coherence of beating for weak pacing
We start by adding stochastic term to Eq. B.11, describing the intrinsic random forces within the cell due to various biological processes: Since we are interested in the long-time coherence of beating, we approximate the short-time noise by a random function χ s (t), which is Gaussian noise with mean χ s (t) = 0 and temporal correlation χ s (t)χ s (t ) ∼ (2D * )δ(t − t )) with D * a measure for the magnitude of the noise. If the external pacing force is weak (A p 1), and the noise varies on a time-scale longer than a single period, we expect that on average the cell will beat with a frequency ω c . To avoid the mathematical complexities of multiplicative noise, we follow the approach described by Hanggi & Riseborough, where a rotating frame of reference is adopted [14]: where y i and y o are the in-phase and out-of-phase components of the beating amplitude respectively. Inserting these expressions into Eq. D.29, and averaging over a single cycle as before yields the two dynamical equations for the amplitude components: withχ s,i (t) andχ s,o (t) the (time averaged) noise in the in and out of phase components with the property χ s,j (t)χ s,j (t ) = 2D * /ω 2 c δ(t − t ) = 2D δ(t − t ), and χ s,j (t)χ s,k (t) = 0 if j = k. The two coupled Langevin equations above are, in general, difficult to solve. We therefore go from the Langevin equations, to the Fokker-Planck (FP) formulation for the probability density P (y i , y o , t; y 0 i , y 0 o , 0), which is the probability that amplitude components y i , y o at time t have certain values, given that at time t = 0 those had the values y 0 i , y 0 o . In vector form, the FP equation is derived [15] from probability conservation, and is given by: with R the n-dimensional drift term vector, corresponding to the non-stochastic part of each of the Langevin equations. Accordingly, for Eq. D.31 we have: We now express these equations in a polar representation where y i = A(t) cos(φ c (t)), y o = A(t) sin(φ c (t)). The FP equation in this representation is: where the drift terms become: Inserting the transformation into Eqs. D.33 and D.35 and re-arranging gives: Where we define: Eq. D.34 still contain terms that couple the amplitude A and phase φ c . To further simplify the problem, we consider the case of weak external force A p 1, and examine the limit of weak noise D /ω c ρ A. For this case, the amplitude has small fluctuations around its steady state value A s = 4ρ/3λω 2 p . Applying these simplifications, and defining D = D /A 2 s and α = α * /A s yields the simplified FP equation described in the main text: with steady-state solution: We further expand Eq. D.39 around its peak φ c = (β + π/2) to get: which for small noise (D → 0) yields: Thus the steady-state distribution for the case of small noise can be approximated by a delta-function of the phase around its steady state value of φ c = β + π/2.

D.2 Dynamics and noise for a strongly paced cell
For a strong external field (pacing), we expect the cell to be entrained to the probe, even in the presence of noise. We therefore use the reference state where the phase φ p (t) is relative to ω p t as in Sec. B (i.e. -χ(t) = A(t) cos(ω p t+φ p (t))). The appropriate FP equation (derived by Hanggi & Riseborough [14] is: For α ∼ A p /A s a measure of the external pacing strength. For a relatively large amplitude of external force, the steady state probability density has been calculated by Hanggi & Riseborough [14]: with Z the normalization and V (φ p ) = − (∆ωφ p − α cos(φ p )) a linearlybiased, periodic potential due to the probe. Eq.D.43 can be further simplified by defining the transformation ψ = φ − φ p . Writing P s in terms of φ p and ψ we have: To proceed, we define the following expansions in terms of the modified Bessel function I n (z) [16]: we calculate the integrals in Eq. D.44,D.46 analytically [16]. This yields the steady-state probability distribution: Fig. 4 we plot Eq. D.47 for various values of the scaled amplitude of probe α/D and detuning ∆ω/D. The probability distribution has the same general characteristics as the stationary probability distribution for the phase of a weakly paced cell presented in Eq. D.39. The distribution is slightly asymmetric about its peak, and this asymmetry is increasing as the detuning ∆ω increases, or as the forcing α decreases. The peak of the distribution shifts towards lower values of φ p as both α and ∆ω increase. For very weak noise, the zero order term in Eq. D.47 is negligible, and one can expand all the modified bessel functions in D 1, and around φ p = π, to yield: which is again the definition of a delta-function δ(φ p − π). The shift of the peak from φ p = (β − π/2) is a result of our shift of the frame of reference from φ c to φ p .
This similarity between the steady-state probability densities allows us to examine the evolution in time of the probability density P (φ, t) by evaluating the moments of the distribution. To this end, we modify the method described by Saito [17], by assuming that the time-dependent probability distribution that satisfies Eq. D.42 has a general form similar to the steady-state distribution of Eq. D.39, but with time dependent coefficients: P (φ, t) = exp (ν(t) cos(φ s − µ(t))) 2π I 0 (ν(t) (D.49) where µ(t) and ν(t) are the two time dependent variables related to peak and width of the probability distribution. By introducing Eq. D.49 into Eq. D.42, and integrating over the phase, we get two dynamical equations for the peak and width:μ = ∆ω − f s cos(µ) I 0 (ν) For small noise, this set of equations yields the steady state expressions: which, for a strongly paced cell (f s ∆ω), is identical to the result in Eq. D.39. This implies that the steady state distribution is mainly determined by the ratio of the reduced amplitude of the external force f s to the noise D.
In the limit of small noise, µ approaches its steady-state value in Eq. D.51 on short time scales of order ∼ 1/f s , and the two dynamical equations in D.50 are decoupled. By expanding For D 1 we can find the approximate expression for the width: This expression is used to estimate the correlation in displacement: x(φ, t)x(φ 0 , 0) = 2π 0 x 0 x(φ, t)P (φ, t) dφ ∼ cos(ω p t+µ s ) 1 − D ν s 1 − e −fst (D.53) where x 0 , φ 0 are arbitrary initial displacement and phase respectively. Note that at long times the correlation function decays to steady state, oscillatory behavior with frequency ω p , and an amplitude smaller than the deterministic case by a term proportional to the noise, D. This implies that as long as the cell is paced by an external force, the phase does not become decorrelated in time, and shows only small fluctuations about its steady-state, oscillatory solution.