Motivated by the need for an analytical tool that can be used routinely to analyze data collected from isolated, detergent-skinned cardiac muscle fibers, we developed a mathematical model for representing the force response to step changes in muscle length (i.e., quick stretch and release). Our proposed model is reasonably simple, consisting of only five parameters representing: (1) the rate constant by which length change–induced distortion of elastic elements is dissipated; (2) the stiffness of the muscle fiber; (3) the amplitude of length-mediated recruitment of stiffness elements; (4) the rate constant by which this length-mediated recruitment takes place; and (5) the magnitude of the nonlinear interaction term by which distortion of elastic elements affects the number of recruited stiffness elements. Fitting this model to a family of force recordings representing responses to eight amplitudes of step length change (±2.0% baseline muscle length in 0.5% increments) enabled four things: (1) reproduction of all the identifiable features seen in a family of force responses to both positive and negative length changes; (2) close fitting of all records from the whole family of these responses with very little residual error; (3) estimation of all five model parameters with a great degree of certainty; and (4) importantly, ready discrimination between cardiac muscle fibers with different contractile regulatory proteins but showing only subtly different contractile function. We recommend this mathematical model as an analytic tool for routine use in studies of cardiac muscle fiber contractile function. Such model-based analysis gives novel insight to the contractile behavior of cardiac muscle fibers, and it is useful for characterizing the mechanistic effects that alterations of cardiac contractile proteins have on cardiac contractile function.

## INTRODUCTION

The recording of the force response of skinned muscle fiber preparations to step-like changes in muscle length (i.e., quick stretch and quick release) has long been a standard experimental procedure in the study of muscle contractile function. Early use of this force response was to interpret the mechanical manifestation of underlying biophysical and biochemical theories of actin–myosin interaction (Huxley and Simmons, 1971; Ford et al., 1977). These and other more recent attempts to represent the force response in terms of mathematical models of muscle based on underlying biochemical events have been reviewed (Kawai and Halvorson, 2007). In addition to validating biophysical theories of muscle contraction, the descriptive features of the force response have proven to be useful for characterizing the contractile function of one muscle as different from another. This is particularly true when this force response has been used for characterizing cardiac muscle (Steiger, 1977; Stelzer et al., 2006a,b, 2007, 2008). Characterization of cardiac muscle behavior based on the force response profile typically involves the fitting of a mono- or bi-exponential function to specific phases of a single step response.

One prominent feature that is missed using single-feature analysis is accounting for the fact that both the overall shape as well as the amplitude of the force response waveform depend on the magnitude of stretch or release, particularly in cardiac muscle. For example, the force responses to large-amplitude stretch and to large-amplitude release are qualitatively very different in shape, whereas the force responses to small-amplitude stretch and to small-amplitude release are similar in shape. This amplitude and directional dependence of the force response demonstrate the existence of a nonlinear contractile feature of cardiac muscle. Such nonlinear features represent an essential aspect of contractile function. Thus, contractile information can be extracted not only by analyzing the force response waveform to a single perturbation, but also by analyzing the overall behavior of the family of force responses to a range of stretches and releases. An appropriately formulated mathematical model could, in principle, capture most of the information contained in the collective family of force responses and allow a general interpretation of the contractile behavior of experimental preparations.

We previously developed a linear mathematical model to describe the force response of constantly activated cardiac muscle to small-amplitude sinusoidal changes in muscle length (Campbell et al., 2004). This model is capable of extracting information pertaining to myofilament contractile dynamic processes and was used as a tool to determine how alterations in cardiac contractile proteins affected myofiber contractile dynamics. For example, we used the model to demonstrate that rat cardiac troponin T (cTnT) modulates sarcomere length-dependent cross-bridge (XB) recruitment (Chandra et al., 2006), and that interaction between myosin heavy chain and troponin isoforms modulates cardiac myofiber contractile dynamics (Chandra et al., 2007).

The linear model consisted of two dynamic variables: (1) a recruitment variable that represents the number of force-bearing XBs; and (2) a distortion variable that represents the average distortion among the population of force-bearing XBs. The recruitment variable varied with muscle length and also varied with the level of Ca^{2+}-dependent activation. The distortion variable underwent transient changes after changes in muscle length. The linear model possessed only four constant-valued parameters: one scaling factor and one rate constant for each of the recruitment and distortion variables. When fit to the data, this simple four-parameter model routinely captured >99% of the variation in force produced by small-amplitude sinusoidal length changes over a broad range of frequencies. All four model parameters were estimated with a high degree of confidence; i.e., there was certainty in the uniqueness of the estimated parameter values. Therefore, this model, when fit to the force response to small-amplitude sinusoidal length perturbations, was recommended as a means for identifying the contractile dynamics of constantly activated cardiac muscle. However, sinusoidal length perturbation analysis is restricted to small-amplitude perturbations and does not elucidate the important nonlinear behaviors described above.

We sought to determine whether our simple linear recruitment–distortion (RD) model could be modified with a nonlinear term for application to characterize the family of force responses to graded amplitude quick stretches and quick releases in constantly activated cardiac muscle. Importantly, could the model be formulated with parameters that allowed effective discrimination between the contractile behaviors of two groups of cardiac muscle bundles whose step responses were not obviously different? We compared step response data from rat cardiac muscle reconstituted with either wild-type rat cTnT (WT-cTnT) or protein kinase C phosphorylation mimetic S199E and T204E mutations in cTnT (cTnT_{S199E/T204E}). Descriptive features of the family of step responses were poorly interpreted and not easily discriminated between the two types of cardiac muscle bundles based on model fits using the simple linear model. Thus, we refined the model to be able to provide meaningful distinctions between these families of force responses.

We found that nonlinearity was an important feature of the step response and that the previous linear model could be expanded with the addition of a nonlinear term to approximate these nonlinear behavioral features. The nonlinear term consisted of a single, one-parameter expression representing interaction between XB distortion and recruitment processes. This single-term addition resulted in a model that could account for most of the nonlinear response. Parameters of this nonlinear RD (NLRD) model were estimated by fitting to data, and these model-estimated parameters provided meaningful distinction between the two groups of cardiac muscle bundles. The remarkable outcome of this work is an NLRD model of constantly activated cardiac muscle that can be used effectively to identify contractile dynamics of muscle bundles responding to step changes in length.

## MATERIALS AND METHODS

### Animal treatment protocols

Heart tissue from which cardiac muscle bundles were dissected was obtained from young adult Sprague-Dawley rats according to a protocol approved by the Washington State University Institutional Animal Care and Use Committee. All animals in this study received humane care in compliance with the animal use principles of the National Academy of Sciences Guide for the Care and Use of Laboratory Animals.

### Preparation of detergent-skinned cardiac muscle fiber bundles

Procedures for the preparation of skinned cardiac muscle fiber bundles have been described elsewhere (Chandra et al., 2007). In brief, rats were anesthetized by inhalation of isoflurane. Hearts were rapidly excised and placed in an ice-cold relaxing solution containing (in mM): 20 BDM, 50 *N*,*N*-bis(2-hydroxyethyl)-2-amino-ethane-sulfonic acid (BES), 20 EGTA, 6.29 MgCl_{2}, 6.09 Na_{2}ATP, 30.83 potassium propionate, 10 sodium azide, 1.0 DTT, and 4 benzamidine-HCl. The pH of the solution was adjusted to 7.0 with KOH. A fresh cocktail of protease inhibitors (in µM: 5 bestatin, 2 E-64, 10 leupeptin, 1 pepstatin, and 200 PMSF) was added to the ice-cold buffered solution. Papillary muscles from these hearts were dissected into bundles ∼0.2-mm wide and 2-mm long. Cardiac muscle bundles were demembranated at 4°C overnight in a relaxing solution containing 1% Triton X-100. Approximately six bundles were studied from each rat heart.

### Reconstitution of recombinant troponin into detergent-skinned rat cardiac muscle fiber bundles

Procedures for reconstitution of endogenous troponin in adult rat cardiac muscle fiber bundles was done based on methods previously described (Chandra et al., 1999, 2006). Endogenous troponin complexes were replaced in skinned cardiac muscle fiber bundles by exposing the fibers to recombinant rat cardiac troponin I (cTnI) and either WT-cTnT or cTnT_{S199E/T204E} variants of rat cTnT for ∼4 h. The extraction solution contained (in mM): 50 BES (pH 7.0 at 20°C), 180 KCl, 10 2,3-butane-dione monoxime, 5 EGTA, 6.27 MgCl_{2}, 1.0 DTT, 0.01% NaN_{3}, and 5 MgATP^{2−}. A fresh cocktail of protease inhibitors was also added. Bundles were rinsed after exposure to recombinant cTnI and cTnT with the extraction solution, and then exposed to recombinant rat cardiac troponin C (cTnC).

### Measurement of force responses to various amplitude quick length changes in detergent-skinned cardiac muscle fiber bundles

Skinned muscle fiber bundles were attached to the arm of a displacement motor (model 308B; Aurora Scientific Inc.) on one end and to a force transducer (AE 801; Sensor One Technologies Corporation) on the other end using aluminum T-clips. Mounted bundles were placed in a 15-µl well, where they were bathed in either activating or relaxing solutions. Sarcomere length, SL, was measured by laser light diffraction methods. Visual display of the first-order laser light diffraction allowed judgment as to the accuracy of the SL assessment and the quality of the preparation. SL in all preparations was set to 2.2 µm before Ca^{2+} activation. Fiber bundle length at SL 2.2 µm was recorded as baseline muscle length (ML). Fiber bundle cross-sectional area was calculated from optical measurements of major and minor axes using an assumption of an elliptical cross section.

Ca^{2+}-activated force was measured at maximal level of activation (pCa 4.3). The amount of free [Ca^{2+}] in solution and corresponding pCa value was determined using methods described previously (Fabiato and Fabiato, 1979). Once steady-state isometric force was achieved, the length of the constantly activated muscle was commanded to change according to the following step length change protocol. A command was given to the motor for the muscle to increase its length by 0.5% of ML in a step-like fashion. This length was maintained for 5 s at which time the muscle was commanded to rapidly return to ML. After another 5 s, a command was given to increase muscle length to 1.0% of ML, which was held for 5 s before rapidly returning to ML. This procedure was repeated for additional step increases in muscle length of 1.5 and 2.0%. All step-like changes in muscle length were essentially completed in 1–2 ms. Measurements of force and muscle length were digitally sampled every 1 ms during the entire procedure. In all cases, force stabilized to a steady-state value within ∼1.5 s after a change in length. The first 1.5 s of the force response to the increase in muscle length was taken as a response to quick stretch, and the first 1.5 s of force response after the return of muscle length to ML was taken as a response to quick release. All force and muscle length records were normalized to their respective values just before length change.

### cDNA clones, expression, and purification of recombinant rat cardiac troponin subunits

Full-length cDNA clones for adult rat cTnI, rat cTnT, and rat cTnC were obtained as described previously (Chandra et al., 2006). The cTnT_{S199E/T204E} gene was created using site-directed mutagenesis techniques. cDNA for WT-cTnT (adult rat) was provided by J.P. Jin (University of Iowa, Iowa City, IA). Rat WT-cTnT, cTnT_{S199E/T204E}, cTnI, and cTnC proteins were expressed in BL21*(DE3) cells and purified as described previously (Chandra et al., 2006).

## RESULTS

Representative families of force responses, *F*(*t*), to four amplitudes of stretch and four amplitudes of release from rat cardiac muscle fiber bundles reconstituted with troponin containing either WT-cTnT or cTnT_{S199E/T204E} are shown in Fig. 1. The force response patterns were qualitatively similar when comparing *F*(*t*) from both groups, which made observable discrimination between the *F*(*t*) from each group difficult. However, careful examination of the profiles of *F*(*t*) from each group reveals subtle differences between the contractile behavior of rat cardiac muscle bundles containing WT-cTnT and those containing mutant cTnT_{S199E/T204E}. Before discussing these differences, we highlight some important features of *F*(*t*) that have been traditionally identified by previous investigators (Ford et al., 1977) and which we observed as well in each *F*(*t*) recorded from both groups.

These features can be illustrated in the force response to large-amplitude quick stretch (Fig. 2). This response exhibited characteristic phases and points that have been previously identified by other investigators (Huxley and Simmons, 1973; Abbott and Steiger, 1977; Ford et al., 1977; Steiger, 1977), with each event interpreted as follows. During steady-state Ca^{2+} activation, muscle bundles reach a steady-state force (*F*_{0}) as a result of increased acto-myosin interactions. Ca^{2+} releases the inhibition of the troponin–tropomyosin complex on acto-myosin interactions, and through subsequent interactions between myofilament proteins, it promotes the formation of several force-bearing XBs. Bound XBs contain elastic regions; power stroke–induced distortion of these elastic regions is the fundamental source of muscle force generation (Huxley, 1957; Huxley and Simmons, 1971; Ford et al., 1977). Rapid stretch of the constantly activated muscle fiber bundle resulted in an immediate increase in force (phase 1), to a value *F*_{1}, due to rapid distortion of the elastic regions of bound XBs. Force then rapidly decays as distorted XBs rapidly detach to reequilibrate into the nondistorted state (phase 2). A slow, but eventual increase in force production (phase 3) then ensues as the increase in muscle length results in a length-induced recruitment of additional XBs into the force-bearing state. Muscle force slightly overshoots (phase 4) and then approaches a steady-state value, *F*_{SS}, as the additional length-recruited XBs equilibrate into the force-bearing state. Some of the characteristics observed in the force response to large-amplitude stretch were also observed in the force response to large-amplitude quick release. Conversely, some of these characteristics of the large-amplitude response were not present or occurred at a different time during the response to release, indicating a departure from linear dependence of the shape of the force response to the magnitude of length change. Thus, there was both linear and nonlinear behavior in the family of force responses of both WT-cTnT– and cTnT_{S199E/T204E}-containing muscle fiber bundles to various amplitude changes in length. We identified and quantified these behaviors as described below.

Linear response features are identified as those whose amplitudes scale linearly with the amplitude of an input. As a consequence, in linear systems, the response to a positive input of given amplitude will be a mirror image of the response to a negative input of the same amplitude. Examples of linear behavior were that of *F*_{1} and *F*_{SS} responses, whose values scaled linearly with the magnitude of the length change (ΔL) (Fig. 3). The slope of the relationship between *F*_{1} and ΔL is an index of the stiffness (E_{D}) of the population of bound XBs at the instance of the length change, and it is proportional to the total number of strongly bound XBs at the time of ΔL. The slope of the relationship between *F*_{SS} and ΔL (E_{R}) is an index of the sensitivity of length-mediated recruitment of additional XBs into the force-bearing state, and it is proportional to the incremental change in the number of force generators in response to ΔL.

Both E_{D} and E_{R} were significantly greater in muscle fiber bundles containing WT-cTnT than in those containing cTnT_{S199E/T204E} (Table I). The increased E_{D} observed in bundles containing WT-cTnT correlates well with an observed increase in *F*_{1} because both E_{D} and *F*_{1} can be correlated to the number of parallel, force-generating XBs bound at steady-state activation before stretch (Campbell et al., 2004).

Nonlinear behavior can be identified as that in which the amplitude or shape of a characteristic does not scale linearly to a given input. Examples of nonlinear contractile behavior were observed in the transient of the force response between *F*_{1} and *F*_{SS}. Such nonlinear behavior is illustrated in the observation that the shape and pattern of the force responses to quick stretch were different than that of the force responses to quick release. For example, the shape of the transient in the force response to large-amplitude quick stretch was different than the shape of the transient in the force response to large-amplitude quick release. The transient to quick stretch possessed a well-defined valley or nadir separating a quick component of the transient from a slow component. In contrast, the transient in a quick-release response did not have a corresponding feature, and the transition from the quick component to a slow component was less well defined (Fig. 2 A). In comparison, the shape of the transient in the force response to small-amplitude quick stretch roughly mirrored that of the transient in the force response to small-amplitude quick release (Fig. 2 B). These observations are a result of the more global observation that the pattern of the shape of the transients in force responses to stretches was different than the pattern of the shape of the transients in force response to releases. For example, the transients in the force responses to stretch each dipped toward one another at a nadir (point *F*_{23} in Fig. 2 A), whereas the shape of the transients in the force responses to releases diverged from another and did not possess a common zenith. Therefore, the shape of the transient between *F*_{1} and *F*_{SS} was considered to have a nonlinear dependence on the magnitude of ΔL because the patterns of these shapes in stretch and release were not mirror images of one another.

Nonlinearities can be further demonstrated by examining how specific periods and phases of *F*(*t*) scale with the magnitude of ΔL. Such analysis was used in a historical study by Ford et al. (1977) to identify and account for nonlinear behavior in skeletal muscle, whereby both the magnitude of tension approached at an early recovery phase (*T*_{2}) and the rate at which tension approached *T*_{2} in the force response to quick release varied with the amplitude of length change. In the present study, we note that the time to reach the dip in responses to stretch, *T*_{23}, and the time to reach 90% of *F*_{SS}, *T*_{90}, in responses to release also varied with the magnitude of ΔL (Fig. 4). This indicates nonlinear behavior in *F*(*t*) because the time course, and thus the shape, of the *F*_{1}-*F*_{SS} transient was different at different magnitudes of ΔL. Both the ΔL-*T*_{23} and the ΔL-*T*_{90} relationships were different in rat cardiac muscle bundles containing WT-cTnT and those containing cTnT_{S199E/T204E}, indicating that the nonlinear behavior was different in these different muscle preparations. Additionally, the magnitude of *F*_{23} scaled curvilinearly with the magnitude of ΔL (Fig. 5). Fig. 5 shows additional distinction between bundles containing the WT-cTnT and those containing cTnT_{S199E/T204E}, further indicating that nonlinear behavior was different in these muscle preparations containing different variants of cTnT.

Characterization of nonlinear behavior, as presented above, only allowed subtle discrimination between the two muscle preparations (note that the scale and magnitude of variation between groups in Figs. 4 and 5 are <10%). Thus, the aim of this study was to characterize the nonlinearities contained in the entire family of force responses to various amplitude step-like length changes in such a way as to easily discriminate contractile behavior of muscle preparations containing different contractile elements. Our goal was to formulate a mathematical model that is able to account for both the linearities and nonlinearities described above, and that is able to discriminate between nonlinear behaviors of different muscle preparations. Both linear and nonlinear features of the family of *F*(*t*) represent important aspects of the contractile behavior of the myofilament. Linear features represent linearly elastic properties related to the number of bound XBs based on ΔL-induced distortion and recruitment of force-bearing XBs, whereas nonlinear features represent other mechanisms of length-mediated contractile activation.

### Model Development

#### Mechanistic considerations.

Model development began with an RD formulation that we used successfully to describe the force response to sinusoidal length change (Campbell et al., 2004; Chandra et al., 2006, 2007). In these sinusoidal force responses, only linear behaviors of the muscle were observed, and a linear combination of independent effects due to distortion and recruitment was used. However, because obvious nonlinear behaviors were observed in the step response, we no longer treated the effects of the recruitment and distortion as independent. Rather, we formulated muscle force to be equal to the product of a recruitment variable, η*(t*), and a distortion variable, *x*(*t*), as is consistent with their physical meaning; i.e., η*(t*) represents the net stiffness of a population of parallel stiffness elements, and *x*(*t*) represents the average elastic distortion among these elements. Thus,

The RD concept arises directly from considering linearly elastic myosin XBs as the elemental force generators in muscle. We refer to earlier publications for the origins of these concepts and for the intellectual bridge between myofilament kinetics and RD fiber bundle dynamics (Campbell et al., 2004; Chandra et al., 2007).

The product of η*(t*) and *x*(*t*) in Eq. 1 to produce force is a nonlinear combination of variables. In our earlier work (Campbell et al., 2004; Chandra et al., 2006, 2007), we linearized Eq. 1 to predict linear small-signal behavior. Here, however, we will start with the product combination of η(*t*) and *x*(*t*), as this nonlinear combination may be important to explaining some of the nonlinear behaviors we observed in the recorded step response. We now examine dynamic expressions for the length-induced changes in η(*t*) and *x*(*t*).

After muscle length change during constant Ca^{2+} activation, net stiffness, η*(t*), changes as more or fewer stiffness elements, XBs, are recruited into or out of the stiffness pool. Because of length-dependent activation kinetics and XB cycling, length-mediated recruitment is a dynamic process requiring time to complete after a change in length. In our earlier work (Campbell et al., 2004; Chandra et al., 2006, 2007), we successfully approximated these recruitment dynamics with the first-order linear differential equation:

where *b* is the recruitment rate constant, β_{0} is a scaling factor, and *l _{d}* is a length at which no recruitment occurs.

Also, muscle length change causes a change in the average elastic distortion, *x*(*t*), of stiffness elements. This change happens immediately; i.e., XB stiffness elements are immediately stretched at the moment the muscle fiber is stretched. However, XB cycling results in the breaking of stretched XBs and the replacement of these broken XBs with new XBs that did not experience the stretch event. We have successfully approximated the dynamics of this transient distortion change with the linear differential equation:

where *c* is the distortion rate constant, *x*_{0} is a baseline isometric distortion imposed on the XB by the power stroke, and

is the time rate of length change or the muscle fiber velocity.

### Qualitative predictions

For an idealized step in which length change occurs instantaneously, analytical solutions of the recruitment and distortion equations can be derived and are helpful in interpreting the response. For an idealized step of magnitude ΔL, the solutions of Eqs. 2 and 3 are, respectively,

These solutions allow calculation of an immediate response as:

and a steady-state response after all transients have died out as:

where both *F*_{1} and *F*_{SS} are normalized with respect to *F*_{0}, *F*_{0} = η_{0}*x*_{0}. From Eq. 6, we see that the starting model predicts that *F*_{1} is linearly dependent on ΔL, even though η*(t*) and *x*(*t*) combine nonlinearly to determine *F*(*t*) in Eq. 1. This model-predicted linear relationship between *F*_{1} and ΔL, which is consistent with what was observed experimentally, demonstrates that distortion by itself is the determinant of the immediate response, and recruitment does not play a role in this very early part of the step response. Further, the starting model also predicts that *F*_{SS} too is linearly dependent on ΔL in spite of the nonlinear character of Eq. 1, which is also consistent with the experimental results. This result demonstrates that recruitment is the determinant of the steady-state response, and stretch-induced distortion does not play a role in the late parts of the response.

The starting model adequately predicts the immediate and steady-state aspects of the step response and now needs to be examined in terms of its ability to predict the transition between these two extremes. Does the nonlinear combination of the recruitment and distortion variables through the multiplication of Eq. 1 produce any of the nonlinear behaviors that were noted in the experimentally observed responses? We examine two cases: (1) the rate constants of recruitment, *b*, and distortion, *c*, are equal; and (2) the rate constant of distortion is 10 times greater than that of recruitment.

The starting model predicts that transition between *F*_{1} of the immediate response and *F*_{SS} of the steady-state response is greatly affected by the relative values of *b* and *c* (see Fig. A1). When the rate constants for distortion and recruitment are equal (*c* = *b*), the transition occurs monotonically with no feature in the trajectory to distinguish various phases of the response. However, when the rate constants differ by an order of magnitude (*c* = 10**b*), the transition trajectory shows two distinct phases: a fast phase as *x*(*t*) quickly recovers to its isometric value and a slow phase as η*(t*) rises (or falls depending on whether stretch or release) to its eventual new steady-state value as a result of length-mediated recruitment effects. Separation of the fast and slow phases produces a nadir in the response to stretch and a zenith in the response to release. Thus, some aspects of the experimentally observed trajectory between *F*_{1} and *F*_{SS} can be recreated with the starting model by separating the values of the rate constants of distortion and recruitment.

However, despite the appearance of a nadir in the trajectory and a clear demarcation of a point in the trajectory that distinguishes the *x*(*t*) recovery from the η*(t*) approach to a new steady state, many of the notable nonlinear features in the experimentally observed step response are absent: (1) for all intents and purposes, the response to a stretch is a mirror image of the response to a release at all step amplitudes, unlike what was observed experimentally where the transition from *F*_{1} to *F*_{SS} was very different for stretches than for releases; (2) trajectories at the various step amplitudes tend to converge on a *F*_{23} transition point in responses to both stretches and releases, whereas transitions in the experimentally observed responses to releases followed widely separated trajectories; (3) although the *F*_{23} transition point occurs earlier in the response when *c* is made progressively larger than *b*, for a given value of *c* and *b* it occurs at the same time in the response, regardless of the step amplitude (not depicted); and (4) the force value at the *F*_{23} transition point scales linearly with step amplitude and never dips below the initial force, no matter how large the difference between rate constants.

In summary, although the starting model generates many overt features of the experimentally observed step response, the nonlinear product combination of η*(t*) and *x*(*t*) in Eq. 1 does not produce any of the prominent or subtle nonlinear features in the experimentally observed step response.

### Model variations to qualitatively reproduce observed nonlinearity

Because the starting model, with its nonlinear combination of η*(t*) and *x*(*t*) in the force calculation (Eq. 1), did not produce any of the nonlinear behaviors that were noted in the experimental results, we examined variants of this starting model in terms of their ability to reproduce aspects of these nonlinear effects. Results from the starting model, especially in the decomposition of the response into its recruitment and distortion components (see Fig. A1), led us to focus on the recruitment equation as the most likely place in the model where variations could bring about improvement. Two general schemes for variation were attempted: (1) variants of Eq. 2 in which velocity had a negative effect on recruitment; and (2) variants in which there was nonlinear interaction between η*(t*) and *x*(*t*).

In brief, all attempts to introduce velocity effects on recruitment were unsuccessful in improving model performance. This was due, in part, to the fact that these velocity effects preferentially and severely degraded the immediate or *F*_{1} response turning the *F*_{1} versus ΔL relationship from a linear shape (which was experimentally observed) into a curvilinear shape (which was not experimentally observed). In contrast, a nonlinear term with an interaction between η*(t*) and *x*(*t*) had a powerful effect on the transition between *F*_{1} and *F*_{SS} without affecting the *F*_{1} and *F*_{SS} aspects of the response, which were already satisfactorily reproduced by the starting model. The challenge became one of identifying the correct form of the η*(t*)–*x*(*t*) interaction term. We approached this challenge as follows.

The interaction term was incorporated into the differential equation for η*(t*) as:

where *f*(*x*) was one of several competing functions of *x*(*t*) (see Table A1). Actually, *x*(*t*) was normalized as

which had the effect of confining distortion effects to situations where *x*(*t*) deviated from its isometric value, *x*_{0} (refer to Appendix for model in normalized form). Further, by dividing this deviation by *x*_{0}, physical units were removed from this term and the interaction coefficient, γ, carried units of s^{−1}, which then allowed γ to be compared with the *b* coefficient. As with the starting model, the equation for distortion in these model variants remained as given by Eq. 3.

Because of the nonlinear interaction term in Eq. 8, an analytical solution for η*(t*) during the transient between *F*_{1} and *F*_{SS} could not readily be obtained. Consequently, we used numerical integration methods to solve for η(*t*) during this transient. However, the steady-state solution for η*(t*) was not affected by the nonlinear term, and the solution for *x*(*t*) was as before (i.e., Eq. 5). Thus, analytical solutions were found for both *F*_{1} and *F*_{SS}, and these were identical to the linear solutions found for the starting model (e.g., Eqs. 6 and 7). As has already been pointed out, these *F*_{1} and *F*_{SS} solutions were consistent with what was observed experimentally.

Model results for three variants of *f*(*x*) are summarized in the Appendix (Table A1). Of these, only one variant was successful in generating the prominent aspects of the experimentally observed nonlinear response:

We demonstrate results from this model variant in the following figures.

The time course of the predicted step response of this model-variant is shown in Fig. 6. The outstanding result is that this model-variant reproduced many of the most prominent nonlinear features of the experimentally observed response transient including:

The shape of the transient of responses to large-amplitude quick stretches was very different than the shape of the transient of responses to large-amplitude quick releases. However, the response to the smallest-amplitude quick stretch roughly resembled the response to the smallest-amplitude quick release.

The trajectory pattern among responses to different ΔL amplitudes was very different for quick stretches than for quick releases, with the trajectories of the four transients to quick stretches at different amplitudes tending to all approach one another at a nadir, while the trajectories of the four transients to quick releases at different amplitudes exhibited no nadirs and remained apart and distinct from one another throughout the transient.

Within the subfamily of the four responses to quick stretches, the value of the force at the time of the nadir in the transient,

*F*_{23}, did tend to converge on a common point and they scaled curvilinearly with the magnitude of length change, and the time to the nadir,*T*_{23}, varied with length change magnitude.Within the responses to quick release, the shape of the response changed from one with an overshoot for the smallest length change to one without an overshoot and a monotonic approach to

*F*_{SS}for the largest length change.

There are only two features of these model predictions that differ from what was seen experimentally: (1) the curvilinearity in the *F*_{23} versus ΔL relationship is in a different direction in the experimental data (curvilinear up with increasing ΔL) than what it is in the model predictions (curvilinear down with increasing ΔL); and (2) there is no overshoot in the response to large-amplitude quick stretches preceding the final approach to *F*_{SS}.

The RD model containing the nonlinear interaction term in the recruitment equation (Eq. 9) reproduced the qualitative nonlinear features observed experimentally in the step response. Thus, it served as the model for validation using the family of force responses to various amplitude quick stretches and quick releases. To make this validation, we fit the model to each of the records of the *F*(*t*) family taken from rat cardiac bundles containing either the WT-cTnT or cTnT_{S199E/T204E} variant. We then compared the fits to the *F*(*t*) family by the simple RD model with those by the NLRD model predictions. Substantial improvement of fits by the NLRD that could be attributed to reproduction of the observed nonlinear traits was taken as evidence of model validity.

### Model validation and application

#### Model fitting and selection methods.

Model predictions were computed using fourth-order Runge-Kutta numerical integration. Computed model predictions were fitted to data records using a Levenberg-Marquardt regression method to minimize sum of square errors (SSE). Results of fitting were obtained not only for the RD and NLRD models described above, but also for all the competing variants of the model shown in the Appendix (Table A2). Comparisons of models were based on several criteria.

#### Criterion 1: Reproduction of the shape of F(t).

An important criterion in model validation was that the model, when fit to data, reproduced the essential shape of *F*(*t*). As shown in Figs. 6 and A1, and Table A1, the starting RD and NLRD models were capable of reproducing the general shape of *F*(*t*). As a qualitative assessment of how well the model was able to reproduce the essential shape, we inspected the likeness of $F\u02c6(t)$ to *F*(*t*) and inspected the time course of the residual errors after data fitting.

#### Criterion 2: Goodness of fit.

Good fit was indicated by low SSE and an R^{2} value close to 1. SSE and R^{2} were determined by the following calculations:

where $F\u02c6(t)$ is the model predicted force response, *F*(*t*) is the measured force response, and $F\xaf$ is the mean value of *F*(*t*). Goodness of fit was also estimated by R^{2} value given by the linear regression of the relationship between *F*(*t*) and $F\u02c6(t)$.

#### Criterion 3: Confidence in the model parameter estimates.

Low standard errors relative to the parameter estimates (coefficient of variation <0.05) suggest that the model parameters were independent in their effects on the predicted force response. A fit resulting in a low standard error of the parameter estimate would suggest that the parameters in the recruitment and distortion equations had independent effects on $F\u02c6(t)$ and were therefore estimated uniquely. Therefore, suitable models were those that, when fit to data, produced low SSE and high R^{2} values and resulted in low standard errors in the estimation of model parameters.

#### Distinguishing between competing models.

From the group of competing models (Tables A1 and A2), the best model was determined based on respective Akaike’s information criterion (AIC). AIC of each model prediction was computed using the following equation:

where *n _{param}* is the number of parameters in the model, and

*n*is the number of observations in the entire

_{obs}*F(t)*data record. AIC was calculated for each model fit and was used to rank the models for each data record. The model that consistently gave the lowest AIC value was considered the best model out of the group of competing models.

#### Model application.

The models were used to describe the contractile behavior of the experimental preparations. Fitted parameters were used to distinguish contractile behavior between rat cardiac muscle bundles containing either WT-cTnT or cTnT_{S199E/T204E}. This was done by comparing parameter estimates obtained in each group using a two-tailed *t* test (α = 0.05).

### Model validation and application results

#### Criterion 1: Reproduction of the shape of F(t).

Representative fits of RD and NLRD models to *F*(*t*) records collected from WT-cTnT– or cTnT_{S199E/T204E}-containing muscle fiber bundles are shown in Fig. 7. The waveforms of the responses in the NLRD predicted responses were qualitatively closer to those in observed data, reproducing the essential phases and nonlinearities of *F*(*t*) described previously, than were those of the RD prediction (Fig. 7). This was substantiated by the observation that the time course of residual errors showed lesser variation when the NLRD model was used to fit data than was observed in the time course of residual errors when the RD model was used to fit data (not depicted). This decrease in residual error due to the addition of the nonlinear interaction term into the recruitment component of $F\u02c6(t)$ suggests that the shape of NLRD predicted $F\u02c6(t)$ was closer to that of *F*(*t*) than was the shape of simple RD model–predicted $F\u02c6(t)$. Therefore, the NLRD model was more accurate in the reproduction of the essential shape of *F*(*t*) than was the RD model, suggesting that the NLRD model is more appropriate for describing data based on criterion 1.

#### Criterion 2: Goodness of fit.

Models were fitted to 11 families of force responses to various amplitude stretches and releases of maximally activated bundles (*n* = 6 from WT-cTnT–containing bundles; *n* = 5 from cTnT_{S199E/T204E}-containing bundles; *n _{obs}* = 8,800 in each

*F*(

*t*)). Both the RD and NLRD models met the first criterion in model selection and validation reasonably well, producing high R

^{2}values and low sum of square residual errors (Fig. 8 and Table II). R

^{2}ranged from 0.881 to 0.966 in simple RD model fits and ranged from 0.971 to 0.993 in NLRD models fits. The NLRD model resulted in lower overall SSE and higher R

^{2}values than did the simple RD model when fit to data records from rat bundles containing either WT-cTnT or cTnT

_{S199E/T204E}variants. The observed decrease in SSE and increase in R

^{2}from the NLRD model fits compared with fits from the RD model was highly significant (Table II). This suggests that the addition of the nonlinear interaction term into the recruitment component of $F\u02c6(t)$ resulted in a better descriptor of the data from both experimental groups because doing so improved the goodness of fit.

The result of the *F(t)* − $F\u02c6(t)$ regressions are shown in Fig. 8. In WT-cTnT–containing muscle fiber bundles, regression fits are $F\u02c6(t)$ = 0.9258*F*(*t*) + 0.0778 for the RD model and $F\u02c6(t)$ = 0.9650*F*(*t*) + 0.0326 for the NLRD model. In cTnT_{S199E/T204E}-containing bundles, regression fits are $F\u02c6(t)$ = 0.8751*F*(*t*) + 0.1291 for the RD model and $F\u02c6(t)$ = 0.9488*F*(*t*) + 0.0493 for the NLRD model. Therefore, the slope of these regressions approached 1 and the intercepts approached 0; slope and intercept were more close to ideal values in the NLRD model than in the RD model (P < 0.0001). Collectively, these results suggest that the addition of the nonlinear interaction term to η*(t*) was appropriate for fitting to data because NLRD resulted in better fits of $F\u02c6(t)$ to *F(t)* and, therefore, satisfied criterion 2 better than the RD model did.

#### Criterion 3: Confidence in the model parameter estimates.

The standard errors of the parameter estimates relative to the respective parameters (coefficients of variation [CoV]) were satisfactorily low, as shown in Table III. For each parameter estimate, the relative standard error was always <5%. This was consistent across fits to all data records using both RD and NLRD models. Because CoV were consistently low for each parameter in each model, we were able to assume that each of the model parameters had independent effects on the overall model-predicted $F\u02c6(t)$. Therefore, both models provided confidence in the uniqueness of the estimated parameter values.

Between the two models, the NLRD model obtained significantly lower CoV of each parameter estimate when compared with those obtained by the RD model, despite the presence of an additional parameter. Over-parameterized models, when fit to data, often result in non-unique solutions and contain parameters that are estimated with little confidence. In this case, however, the addition of the nonlinear interaction term was not damaging to our confidence in parameter estimation, but instead increased the confidence by which parameters were estimated.

#### Discriminating between competing models.

Because both models to a large degree satisfied criteria 1–3, we looked at the AIC to provide an objective index to rank the competing models (Tables A1 and A2). In every fit of *F*(*t*) from both WT-cTnT– and cTnT_{S199E/T204E}-containing cardiac muscle fiber bundles, the AIC index was less for the NRLD model than it was for the RD model. Based on this index, the NLRD model was always a better descriptor of *F*(*t*) than was the RD model.

#### Model application.

Fitted model parameter estimates from RD and NLRD are shown in Figs. 9 and 10, respectively. The RD model was unable to distinguish any difference in the rate constant of length-mediated XB recruitment, *b*, between groups containing either WT-cTnT or cTnT_{S199E/T204E} (Fig. 9 A). The RD model was, however, capable of distinguishing a faster rate constant of strain-induced XB detachment, *c*, in cardiac muscle fiber bundles containing the cTnT_{S199E/T204E} when compared with muscle fiber bundles containing WT-cTnT (Fig. 9 B).

Similarly, the NLRD model was not able to distinguish differences in the XB recruitment rate constant between muscle fiber bundles containing either WT-cTnT or cTnT_{S199E/T204E}, but it readily distinguished a faster dynamic rate constant of strain-induced XB detachment in muscle fiber bundles containing cTnT_{S199E/T204E} (Fig. 10, A and B). However, the NLRD model provided further distinction between the contractile behavior of the two muscle groups in that the NLRD model predicted that XB recruitment was affected differently in the two groups of muscle fibers by the XB strain-mediated effect on the recruitment of strong XBs (i.e., γ; Fig. 10 C). XB strain had a much greater negative impact on XB recruitment in bundles containing the cTnT_{S199E/T204E} than in muscle fiber bundles containing the WT-cTnT, contributing to a more pronounced nonlinearity in *F*(*t*) that was observed in bundles containing cTnT_{S199E/T204E}. Although this unique characteristic of the cTnT_{S199E/T204E}-containing group was easily identifiable using the NLRD model, the RD model was unable to distinguish between the nonlinear behaviors of groups containing either WT-cTnT or cTnT_{S199E/T204E}.

Therefore, the NLRD model can be used to provide a more detailed interpretation of the effects that alterations in the structure of contractile proteins have on the contractile function of cardiac muscle fiber bundles. The NLRD model contains the nonlinear interaction term, which is formulated as an effect by which the state of bound XBs influences the recruitment of other XBs. Because XBs do not directly interact with one another, this interaction describes cooperative and/or allosteric mechanisms translated along the thick or thin filaments. Here, our study shows that the thin filament is involved in the transduction of the XB strain-dependent effect on other XBs because different variants of cTnT affect this nonlinear interaction process differently. Modifying cTnT and, in turn, the thin filament structure modified the mechanisms by which XBs influence the recruitment of other XBs, as measured by γ. This suggests that thin filament proteins (e.g., cTnT) play an important allosteric regulatory role in XB-mediated XB recruitment, a process that underlies prominent nonlinear contractile behavior in cardiac muscle.

## DISCUSSION

Motivated by the need for an analytical tool that can be used routinely to analyze data collected from isolated, detergent-skinned cardiac muscle fibers, we developed a mathematical model for representing the force response to step changes in muscle length (i.e., quick stretch and release). Fitting this model to a family of force recordings representing responses to eight amplitudes of step length change (±2.0% ML in 0.5% increments) enabled unambiguous estimation of five characteristic contractile parameters that were used to distinguish muscle fibers exhibiting even subtly different contractile characteristics. These model-based differences provide a dynamical contraction profile of the muscle fiber that goes beyond what can be obtained from the suite of information collected from other routine assays in this preparation, such as force-pCa curves, *k*_{tr}, and ATPase-pCa curves. Collectively, with information from the other routine assays, this model-based analysis gives unprecedented insight to the contractile status of the muscle fiber.

The model proposed here differs from an earlier model that successfully represented the force response to small-amplitude sinusoidal length change (Campbell et al., 2004; Chandra et al., 2007) by the addition to that earlier model of a single nonlinear term. This nonlinear term is capable of reproducing nonlinearity observed in the shape of force responses to various amplitude quick stretches and quick releases. This nonlinearity became more evident in large-amplitude step responses but was less evident in small-amplitude step and sinusoidal responses.

The current nonlinear model (and the earlier linear model) was constructed using a combination of mechanistic notions and phenomenologic methods. Mechanism dictated the general structure of the model, which came from the understanding that contractile force depends on the number of parallel stiffness elements and their elastic distortion. Rapid changes in muscle length elastically distort existent stiffness elements to immediately increase the force (the incident phase of the response). These distorted elements are short-lived as they dissipate and are replaced by newly formed elements that did not experience the distortion of rapid length change. Thus, the incident phase due to elastic distortion dissipates quickly, and this imparts a rapid recovery phase to the force transient in the direction of the initial force before length change. However, changes in muscle length recruit more stiffness elements (or fewer, depending on the direction of length change) into the stiffness-generating pool of elements that are responsible for contractile force. This recruitment of stiffness elements occurs more slowly than the dissipation of distortion and imparts a relatively slow dynamic to the transient in the direction away from the initial force. Thus, the model’s general mathematical structure should account for the dynamics of distortion and recruitment during the response.

Phenomenologic methods, which are data driven, determine model size and specific structural form. Thus, the essence of both the incident and the rapid recovery phases of the response appeared as if it could be captured by a single, linear, first-order differential equation whose driving input was the velocity of length change. Such phenomenology was consistent with the mechanistic idea that elastic distortion of replaceable elements was responsible for these phases. In a similar manner, the essence of the slower phase of the response, as it approached its eventual steady state, appeared as if it could be reproduced by a single, linear, first-order differential equation whose driving input was the length change. This was also consistent with the mechanistic idea that recruitment of stiffness elements was responsible for this phase of the response.

In the linear RD model, distortion as an outcome of the distortion equation and recruitment as an outcome of the recruitment equation were independent of one another. However, such independence of these variables was not consistent with the observed nonlinearity indicated by the fact that the time course of the response to large-amplitude stretch was markedly different than the time course of the response to large-amplitude release.

There are many ways to introduce nonlinearity into a set of equations with two variables. Perhaps the simplest is to add a term to the equation for one of the variables in which there is multiplicative interaction between both variables. We introduced such a term into the differential equation for recruitment by negatively coupling recruitment to the square of distortion, making this effect independent of the sign on distortion. The effect of this nonlinear term was to bring about a transient depression in the pool of stiffness elements during the time of the distortion transient whether the distortion was positive or negative. The resulting nonlinear model satisfactorily reproduced the essence of the observed difference between large-amplitude stretches and releases.

Our efforts to build a mathematical model to represent the force response of constantly activated muscle to quick length change differ from other similar efforts (for review see Kawai and Halvorson, 2007) on two accounts. One, our goal was to build an easily applicable analytic tool for routine data analysis of routine experiments. This differs with the goals of others who have sought to build a mathematical bridge between the many biochemical/biophysical steps in the XB cycle and the force response. Two, we used only general notions about underlying mechanisms to inform our model and relied on reproductions of observed phenomenology in the force response to dictate the specific form and content of the final model. This differs with the mechanistic methods in which mathematical models are constructed to account for underlying assumptions that are made with regard to kinetic steps in the XB cycle and to incorporate any number of possible influences of these kinetic steps. Models developed from such mechanistic methods yield valuable insights to the biophysics of muscle contraction but, due to their nature, are often over parameterized and not easily applied in laboratory settings where easily applicable analytical tools are required to interpret experimentally obtained data.

In contrast, our proposed NLRD model is reasonably simple, consisting of only five parameters representing: (1) the rate constant by which length change–induced distortion of elastic elements is dissipated; (2) the stiffness of the muscle fiber; (3) the amplitude of length-mediated recruitment of stiffness elements; (4) the rate constant by which this length-mediated recruitment takes place; and (5) the magnitude of the nonlinear interaction term by which distortion of elastic elements affects the number of recruited stiffness elements. This model reproduces all the identifiable features seen in a family of force responses to both positive and negative length changes. It fits the records of the whole family of these responses with very little residual error, and the parameters of the model are identified with great certainty. Finally, yet quite importantly, application of the model to data recorded from cardiac muscle fibers with different contractile regulatory proteins allowed easy discrimination between muscle fibers with altered contractile function.

For instance, the model distinguished that muscle fiber bundles containing the cTnT_{S199E/T204E} variant exhibited faster XB detachment dynamics (higher *c*) and a greater effect of strong XB strain on the recruitment of other strong XBs (higher γ) when compared with fibers containing WT-cTnT. These model-driven observations can be used for providing mechanistic explanations of physiological observations that may otherwise be unclear. For example, changes in the dynamic behavior of XB recruitment and detachment events may explain the decrease in maximum force generation (Fig. 1, legend) and stiffness (Fig. 3 and Table I) of muscle fibers containing the cTnT_{S199E/T204E} variant. Faster XB detachment dynamics (*c*) and a greater negative effect of strained XBs on the state of other bound XBs (γ) indicate that strong XBs more readily dissociate from thin filaments containing cTnT_{S199E/T204E}. This enhanced dissociation of strongly bound XBs would shift the equilibrium of XB recruitment away from the force-bearing pool, and it gives rise to a lower capacity of force generation and contributes to a lower stiffness in muscle fiber bundles containing cTnT_{S199E/T204E}. Because cTnT_{S199E/T204E} is a phosphorylation analogue of cTnT (by protein kinase C), these model-driven predictions using experimentally obtained data suggest a novel mechanistic explanation for the down-regulation of contractile function by protein kinase C–mediated phosphorylation of cTnT observed by us and by others (Burkart et al., 2003; Sumandea et al., 2003).

Thus, the nonlinear recruitment distortion model is useful for cardiac muscle experimentalists to characterize and understand the effects that alterations in contractile proteins have on cardiac contractile function. This model, when combined with physiological studies, provides mechanistic insights with a high level of confidence to describe contractile function beyond that achieved by most current methods of analysis. Therefore, we recommend this mathematical model as an easily applicable analytic tool for routine use in studies of cardiac muscle fiber contractile function.

## APPENDIX

### Model normalization

The basic RD model takes force to be the product:

where

Re-expressing the RD model in terms of force and length normalized to initial force, *F*_{0}, and initial length, *l*_{0}, gives for an initial force, *F*_{0} = η_{0}*x*_{0}. Then, normalized force is given by

And normalized fiber length is

These allow definition of a normalized stiffness:

And a normalized distortion:

Further definitions can be made of normalized parameters:

Note that from Eq. A2, $\eta 0=[\beta b](l0\u2212ld)$, and thus,

Restating the model in terms of normalized variables:

where the normalized differential equation for XB stiffness becomes:

This can be expressed in terms of normalized muscle length by substituting *l*(*t*) = λ(*t*)*l*_{0} and *l*_{d} = λ_{d}*l*_{0} into the above equation, as follows:

Similarly, the normalized differential equation for the XB distortion component becomes:

Substituting $dl(t)dt=l0d\lambda (t)dt$ (because $ddt\lambda (t)=ddtl(t)l0$) into the above equation,

where the only terms carrying physical units are *b* and *c* with units of s^{−1}.

In the normalized model, there are four parameters to estimate: *b*, λ* _{d}*,

*c*, and υ

*. However, parameter λ*

_{0}*can be estimated from steady-state ϕ*

_{d}_{ss}versus Δλ relationship, as described below. Therefore, only three parameters need to be estimated from fitting to normalized step response transients:

*b*,

*c*, and υ

_{0}.

Starting values for variables in numerical integration: ϕ_{0} = 1, ε_{0} = 1, and ζ_{0} = 1.

### Advantages of expressing in terms of normalized force:

It eliminates the parameter, β, from equations and thus reduces by one the number of unknown parameters.

- It simplifies the estimation of the parameter, λ
, which can easily be done from steady-state ϕ_{d}_{ss}versus Δλ regression line. This relationship is shown at times long after the instance of stretch, where the transient phases of the force response die out and only the length-mediated recruitment component remains:$\phi ss=1+\Delta \lambda 1\u2212\lambda d$$(1\u2212\lambda d)(\phi ss\u22121)=\Delta \lambda $This further reduces the number of free parameters to be estimated from dynamic transients from 4 to 3.$\lambda d=1\u2212\Delta \lambda \phi ss\u22121=\Delta \lambda +1\u2212\phi ss1\u2212\phi ss$ - It simplifies specifying initial value for υ
_{0}estimation as follows. Let ϕ_{1}be the peak of immediate response. Then, to an approximation,and a close approximation of υ$\phi 1=1+\Delta \lambda \upsilon 0$_{0}can be obtained from the slope of the ϕ_{1}versus Δλ regression line. This approximation can then be used as a starting value for the υ_{0}estimation. It simplifies specification of starting values for variables, ε

*(t*) and ζ*(t*), when integrating numerically.

### Step response predictions of normalized RD model

#### Accounting for residuals.

Assuming that the residuals have a linear dependence on ΔL, we can account for the non-normal distribution of the residuals by modifying the model as follows:

As was the case with the other models variants, this model variation improved fitting by further minimizing residual SSE, but did not provide strong statistical (i.e., AIC) measures to warrant its use over the NLRD model.

## Acknowledgments

This work was supported, in part, by the National Heart, Lung, and Blood Institute (grant R01-HL75643 to M. Chandra and grant R01- HL80186 to W. Dong), ARCS Fellowship (to S.J. Ford), and Poncin Fellowship (to R. Mamidi).

Richard L. Moss served as editor.