Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.12(6); 2016 June**|**PMC4898704

Formats

Article sections

Authors

Related links

PLoS Comput Biol. 2016 June; 12(6): e1004904.

Published online 2016 June 8. doi: 10.1371/journal.pcbi.1004904

PMCID: PMC4898704

Andrew D. McCulloch, Editor^{}

University of California San Diego, UNITED STATES

The authors have declared that no competing interests exist.

Conceived and designed the experiments: TL WH. Performed the experiments: TL. Analyzed the data: GST TL GD WH. Contributed reagents/materials/analysis tools: GST TL. Wrote the paper: GST WH.

* E-mail: ta.zarg-inu@plit-rehcappahcs.nurdug

Received 2015 October 1; Accepted 2016 April 5.

Copyright © 2016 Schappacher-Tilp et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Passive forces in sarcomeres are mainly related to the giant protein titin. Titin’s extensible region consists of spring-like elements acting in series. In skeletal muscles these elements are the PEVK segment, two distinct immunoglobulin (Ig) domain regions (proximal and distal), and a N2A portion. While distal Ig domains are thought to form inextensible end filaments in intact sarcomeres, proximal Ig domains unfold in a force- and time-dependent manner. In length-ramp experiments of single titin strands, sequential unfolding of Ig domains leads to a typical saw-tooth pattern in force-elongation curves which can be simulated by Monte Carlo simulations. In sarcomeres, where more than a thousand titin strands are arranged in parallel, numerous Monte Carlo simulations are required to estimate the resultant force of all titin filaments based on the non-uniform titin elongations. To simplify calculations, the stochastic model of passive forces is often replaced by linear or non-linear deterministic and phenomenological functions. However, new theories of muscle contraction are based on the hypothesized binding of titin to the actin filament upon activation, and thereby on a prominent role of the structural properties of titin. Therefore, these theories necessitate a detailed analysis of titin forces in length-ramp experiments. In our study we present a simple and efficient alternative to Monte Carlo simulations. Based on a structural titin model, we calculate the exact probability distributions of unfolded Ig domains under length-ramp conditions needed for rigorous analysis of expected forces, distribution of unfolding forces, etc. Due to the generality of our model, the approach is applicable to a wide range of stochastic protein unfolding problems.

We provide a simple and stable algorithm to determine the exact solution of passive forces in a half sarcomere in length-ramp simulations. The approach is applicable to a wide range of stochastic models of protein unfolding.

Passive forces in sarcomeres or myofibrils are almost exclusively governed by the giant protein titin [1]. A titin strand spans the half sarcomere from Z-disk to M-band. While its section located in the thick filament is nearly inextensible, its I-band region functions as a molecular spring. In skeletal muscles, titin’s I-band region comprises two immunoglobulin (Ig) domains, a N2A portion as well as a region rich in proline(P), glutamate (E), valine (V), and lysine(K), the PEVK region [1, 2]. The distal Ig domains (close to the AI-junction) are thought to form almost inextensible end-filaments [3, 4] (Fig 1).

The proximal Ig domains (close to the Z-disk) are able to unfold in a force and time dependent manner [5]. Single molecule atomic force microscopy, in which the length of a single molecule is controlled, show a characteristic sawtooth pattern in force-extension curves due to a complex time series of unfolding events, e.g. [6–10]. Molecular dynamic simulations coincide with experimental data and provide new insight into the mechanics of unfolding at the atomic level; e.g. [11–14].

In a half sarcomere around 1.210^{9} titin strands per mm² are arranged in parallel [15, 16]. Therefore, the sawtooth pattern in half-sarcomers is completely averaged out. However, it is still possible to observe unfolding events in myofibrils implicitly because of a prominent change of stiffness [17, 18].

A favourable experimental set-up to study unfolding characteristics, like mean dwell times or the distribution of forces at rupture of single proteins, are so-called force-clamp studies and force-ramp studies. In force-clamp studies the force is maintained at a constant level whereas in force-ramp studies the force increases linearly [19]. Elongation of protein length due to unfolding events is then a step function over time [20, 21]. Experimental data can be rigorously analysed by order statistics [22, 23].

However, the theoretical framework of force-clamp analysis is not applicable in length-ramp experiments which are frequently used in myofibril studies, like [24–28] to name a few. Models of active and passive force production under dynamic conditions on the sarcomere or myofibril level either use Monte Carlo simulations of single titin strands where single protein forces are scaled to half sarcomere forces, e.g. [29], or phenomenological models, e.g. [30, 31]. Monte Carlo simulations of single titin strands deliver detailed information about the possible behaviour of single titin strands which is not necessarily needed for models on the half sarcomere level. On the other hand the characteristics of force production on the sarcomere level are still only an approximation. Phenomenological models focus on the characteristics of passive forces in half sarcomeres thereby condensing the information but they might not be sufficient when dealing with new theories on muscle contraction involving titin binding to actin upon activation (e.g. [25, 29, 32, 33]). These theories require knowledge of titin’s mechanical behavior as well as the mechanics of isolated titin segments

In this work we are interested in the forces contributed by all titin strands in a half sarcomere. Since titin strands are aligned in parallel the total force is the sum of single protein forces and thereby the weighted expectation value of a single molecule. We reformulate a well known titin model which is based on the kinetic properties of titin’s macroscopic structures and therefore easily adjustable to various theories of muscle contraction. We define a simple two state unfolding regime and derive the exact probabilities describing the regime. Finally, we compare our results to Monte Carlo approximations and previously published results, as well as our own experimental data. Our approach provides a simple algorithm to calculate the exact expectation value of forces contributed by titin in a half sarcomere or probability distribution of unfolding forces.

Throughout the manuscript we base our simulations on rabbit psoas 3.400-kD isoform [34], i.e. 50 proximal Ig domains, 800 PEVK residues and 26 distal Ig domains. However, the isoform only affects numerical results but not the theoretical considerations.

The stochastic model in a slightly different form, and the corresponding Monte Carlo simulations have been published before, e.g. [7, 29, 35–37]. In this work, we provide a rigorous but simple mathematical formulation of the model allowing for a straightforward deduction of the existence of a unique solution of the stochastic model. In addition, we include the possible hierarchy of unfolding [36] or (mathematically equivalent) different mechanical stability of Ig domains [38] in both, Monte Carlo approximations and the exact solution.

To allow an effective description of the theoretical considerations, we formulate the theoretical framework for stretching experiments where refolding is highly unlikely. Finally, we introduce refolding events for a simplified case.

Single-molecule experiments suggest that titin’s Ig domains (folded or unfolded) act like a molecular spring showing wormlike chain (WLC) behaviour [9]. Upon unfolding the extended Ig domain gains the distance *d*_{u}. Titin’s PEVK region shows modified WLC behavior [16]. A well known approximation to the WLC model of entropic elasticity [39, 40] relates the external force *F* to the end to end length *x* of the chain by

$$\begin{array}{c}\hfill F=\frac{{k}_{B}T}{pl}\left(\frac{1}{4{(1-\frac{x}{cl})}^{2}}-\frac{1}{4}+\frac{x}{cl}\right),\end{array}$$

(1)

where *pl* is the persistence length and *cl* the chain’s contour length, *k*_{B} is the Boltzmann constant and *T* the absolute temperature. The modified WLC model includes an enthalpic contribution to elasticity by including a stretch modulus *F*_{0}:

$$\begin{array}{c}\hfill F=\frac{{k}_{B}T}{pl}\left(\frac{1}{4{(1-\frac{x}{cl}+\frac{F}{{F}_{0}})}^{2}}-\frac{1}{4}+\frac{x}{cl}-\frac{F}{{F}_{0}}\right).\end{array}$$

(2)

The stochastic model describing the force contributed by a single titin strand in a half sarcomere with length *l* can be expressed as a convex optimization problem:

$$\begin{array}{c}\hfill \begin{array}{c}\underset{{l}_{P\phantom{\rule{-0.166667em}{0ex}}E\phantom{\rule{-0.166667em}{0ex}}V\phantom{\rule{-0.166667em}{0ex}}K},{l}_{f},{l}_{u},{l}_{e\phantom{\rule{-0.166667em}{0ex}}n\phantom{\rule{-0.166667em}{0ex}}d},{l}_{r\phantom{\rule{-0.166667em}{0ex}}e\phantom{\rule{-0.166667em}{0ex}}s\phantom{\rule{-0.166667em}{0ex}}t}}{\text{minimize}}\left({V}_{f}\left({l}_{f}\right)+{V}_{u}\left({l}_{u}\right)+{V}_{P\phantom{\rule{-0.166667em}{0ex}}E\phantom{\rule{-0.166667em}{0ex}}V\phantom{\rule{-0.166667em}{0ex}}K}\left({l}_{P\phantom{\rule{-0.166667em}{0ex}}E\phantom{\rule{-0.166667em}{0ex}}V\phantom{\rule{-0.166667em}{0ex}}K}\right)+{V}_{end}\left({l}_{e\phantom{\rule{-0.166667em}{0ex}}n\phantom{\rule{-0.166667em}{0ex}}d}\right)+{V}_{r\phantom{\rule{-0.166667em}{0ex}}e\phantom{\rule{-0.166667em}{0ex}}s\phantom{\rule{-0.166667em}{0ex}}t}\left({l}_{r\phantom{\rule{-0.166667em}{0ex}}e\phantom{\rule{-0.166667em}{0ex}}s\phantom{\rule{-0.166667em}{0ex}}t}\right)\right)\hfill \\ \text{subject to}\left\{\begin{array}{c}l={l}_{P\phantom{\rule{-0.166667em}{0ex}}E\phantom{\rule{-0.166667em}{0ex}}V\phantom{\rule{-0.166667em}{0ex}}K}+{l}_{f}{n}_{f}+{l}_{u}(N-{n}_{f})+{l}_{e\phantom{\rule{-0.166667em}{0ex}}n\phantom{\rule{-0.166667em}{0ex}}d}+{l}_{r\phantom{\rule{-0.166667em}{0ex}}e\phantom{\rule{-0.166667em}{0ex}}s\phantom{\rule{-0.166667em}{0ex}}t}\hfill \\ {l}_{P\phantom{\rule{-0.166667em}{0ex}}E\phantom{\rule{-0.166667em}{0ex}}V\phantom{\rule{-0.166667em}{0ex}}K},{l}_{f},{l}_{u},{l}_{e\phantom{\rule{-0.166667em}{0ex}}n\phantom{\rule{-0.166667em}{0ex}}d},{l}_{r\phantom{\rule{-0.166667em}{0ex}}e\phantom{\rule{-0.166667em}{0ex}}s\phantom{\rule{-0.166667em}{0ex}}t}\ge 0\hfill \end{array}\right.,\hfill \end{array}\end{array}$$

(3)

where *l*_{PEVK}, *l*_{f}, *l*_{u}, *l*_{end} and *l*_{rest} are the lengths of PEVK, folded and unfolded Ig domains, end-filaments as well as the sum of titin’s A-band region and slack length respectively. *V*^{u, f} is the potential of the WLC model for unfolded and folded proximal Ig domains, *V*_{PEVK} is the potential of the modified WLC model. *V*_{end} and *V*_{rest} are potentials of very stiff linear springs describing the almost in-extensible end-filaments formed by distal Ig domains, and titin’s A-band region and slack length, respectively. *N* is the total number of Ig domains. The number of folded Ig domains *n*_{f} is a discrete random variable. The formulation of the optimization problem guarantees the existence of a unique solution of Eq (3) whenever the number of folded Ig domains is known and the convex feasible set is non-empty.

Mechanical model parameters are average values obtained from the literature and provided in Table 1.

Since we are interested in simulating length-ramp experiments, we solve Eq (3) for every possible *n*_{f} (see Fig 2) for half sarcomere lengths starting at 1μm. With *F*_{i}(*l*(*t*)) = *F*_{i}(*t*), *i* = 0, …, 50 we denote the force as a function of length (and hence time) when exactly *i* Ig domains are unfolded.

Unfolding and refolding of one Ig domain is a thermally driven process described by the unfolding rate function *u* and the folding rate function *f*. The barrier height of the unfolding and refolding process is influenced by the external force *F*:

$$\begin{array}{c}\hfill u\left(F\right)={\omega}_{0}{e}^{\frac{F{x}_{u}}{{k}_{b}T}},r\left(F\right)={\omega}_{1}{e}^{-\frac{F{x}_{f}}{{k}_{b}T}},\end{array}$$

(4)

where *ω*_{0} and *ω*_{1} are spontaneous unfolding / folding rates at zero force and *x*_{u}, *x*_{f} are the width of the activation barriers, e.g. [35]. Due to the length-ramp condition the force changes significantly with every unfolding event thereby directly influencing the unfolding probability of the remaining folded Ig domains.

Since there is experimental support for hierarchical unfolding of Ig domains [36], or differences in the mechanical stability of Ig domains homogeneously distributed across the titin strands [38], we introduce *m* Ig clusters, where each cluster has a different spontaneous unfolding rate *ω*_{0}. We denote the different unfolding rates by *u*_{k}, *k* = 1, …, *m* and assume (without loss of generality) that *u*_{1}(*F*) ≤ *u*_{2}(*F*) ≤ … ≤ *u*_{m}(*F*) for any given force *F*. It is worth pointing out that we could define the clusters by different activation barriers if it would prove favourable in terms of parameter estimation or data fitting. However, such an approach would not change the characteristics of the simulations presented in this paper.

We simulate *R* titin strands and approximate the expected forces in a half sarcomere by averaging the forces of *R* titin strands and scaling the force by the number of titin strands in a half sarcomere. For a single titin strand let *m* be the number of unfolding clusters. We simulate stretches starting at low half sarcomere lengths where unfolding events are highly unlikely. Therefore, we start with 50 folded Ig domains. The probability *p*_{k} that one Ig domain of the *k*-th cluster with *n*_{k} Ig domains will unfold within the next time step is approximately *p*_{k}(*l*(*t*_{1})) = *n*_{k} *u*_{k}(*F*_{0}(*l*(*t*_{1}))) *dt*. By drawing an equally distributed random number 0 ≤ *z* ≤ 1, we define an unfolding event in the *k*-th cluster if

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \sum _{j=1}^{k-1}{p}_{j}\left(l\left({t}_{1}\right)\right)& <z\le & {p}_{k}\left(l\left({t}_{1}\right)\right)\hfill \\ \hfill \text{or\hspace{0.17em}}\sum _{j=1}^{k-1}{n}_{j}\xb7{u}_{j}\left({F}_{0}\left({t}_{1}\right)\right)\xb7dt& <z\le & {n}_{k}\xb7{u}_{k}\left({F}_{0}\left({t}_{1}\right)\right)\xb7dt\hfill \end{array}\end{array}$$

(5)

If no cluster fulfils condition (5), no unfolding event takes place. In a similar way, we determine an unfolding event at any other time step *t*_{i}, where *k*_{1}, *k*_{2}, …, *k*_{m} Ig proteins of the first, second,…,*m*th cluster are already unfolded, respectively. Again, let 0 ≤ *z* ≤ 1 be an equally distributed random number. Let ${N}_{f}=N-{\sum}_{j=1}^{m}{k}_{j}$. We define an unfolding event in the *l*-th cluster if

$$\begin{array}{c}\hfill \sum _{j=1}^{l-1}({n}_{j}-{k}_{j})\xb7{u}_{j}\left({F}_{{N}_{f}}\left({t}_{i}\right)\right)\xb7dt<z\le ({n}_{l}-{k}_{l})\xb7{u}_{l}\left({F}_{{N}_{f}}\left({t}_{i}\right)\right)\xb7dt\end{array}$$

(6)

If refolding has to be taken into account, e.g., simulation of hysteresis curves [5, 17, 18], the formula can be adjusted accordingly. For the sake of simplicity, we introduce the formula for one cluster which can be easily generalized to an arbitrary amount of clusters: Let *i* Ig domains be unfolded at some time *t*. The probability that a folded Ig domain will unfold within the next time step *t* + *dt* is approximately (*n* − *i*) *u*(*F*_{i}(*l*(*t* + *dt*))) *dt*, while the probability that one of the *i* unfolded Ig domains will refold within the next time step is approximately *i* *r*(*F*_{i}(*l*(*t* + *dt*))) *dt*. By drawing an equally distributed random number 0 ≤ *z* ≤ 1, we define an unfolding event if

(7)

and a refolding event if

(*n* - *i*) · *u*(*F*_{i}(*l*(*t* + *d**t*))) · *d**t* < *z* ≤ (*n* - *i*) · *u*(*F*_{i}(*l*(*t* + *d**t*))) · *d**t* + *i* · *r*(*F*_{i}(*l*(*t* + *d**t*))) · *d**t*.

(8)

Finally, we determined the number *R* of simulated titin strands by estimating the error of Monte Carlo simulations. Specifically, we performed *R* independent Monte Carlo simulations *H* times to estimate the mean force *F*_{h}(*l*_{i}) at a given length *l*_{i}, *i* = 1, …, *L* for *h* = 1, …, *H*. We define the (*R* dependent) error of the Monte Carlo simulation *M**C**E* by

$$\begin{array}{c}M\text{\hspace{0.17em}}C\text{\hspace{0.17em}}E{\left(R\right)}^{2}={\Vert \langle {F}_{h}^{2}(\xb7)\rangle -{\langle {F}_{h}(\xb7)\rangle}^{2}\Vert}_{\infty},\end{array}$$

where $\langle x\rangle =(1/N)\xb7{\sum}_{n=1}^{N}{x}_{i}$ for *x* = (*x*_{1}, …, *x*_{n}). For *H* = 100 we get the following *M**C**E* estimates: *M**C**E*(10) = 2.71, *M**C**E*(50) = 1.3, *M**C**E*(100) = 0.84, *M**C**E*(200) = 0.63, *M**C**E*(500) = 0.39. Therefore, seeking a compromise between simulation time and accuracy, we chose to simulate 200 titin strands.

For the exact solution, we need to calculate the probabilities that after a stretch to a given length *l*(*t*) no unfolding event took place (*P*_{0}), exactly one unfolding event took place (*P*_{1}), up to *N* unfolding events took place (*P*_{N}). The mathematical formulation of the probabilities is straight forward. For instance, *P*_{0} is the survival probability of all Ig domains, i.e.

$$\begin{array}{c}\hfill {P}_{0}\left(t\right)=exp\left(-\sum _{k=1}^{m}{n}_{k}{\int}_{0}^{t}{u}_{k}\left({F}_{0}\left(\tau \right)\right)d\tau \right).\end{array}$$

(9)

*P*_{1}(*t*) is the probability, that no unfolding takes place until time *τ*, then one Ig domain of one of the *m* clusters unfolds, and then no further event takes place between *τ* and *t*, i.e.

$$\begin{array}{c}\hfill {P}_{1}\left(t\right)=\sum _{l=1}^{m}{\int}_{0}^{t}{P}_{0}\left(\tau \right){u}_{k}\left({F}_{0}\left(\tau \right)\right)exp\left(-\sum _{i=1}^{m}({n}_{i}-{\delta}_{i,l}){\int}_{\tau}^{t}{u}_{i}\left({F}_{1}\left(\sigma \right)\right)d\sigma \right)d\tau ,\end{array}$$

(10)

where ${\delta}_{ij}=\{\begin{array}{cc}\hfill 0,& i\ne j\hfill \\ \hfill 1,& i=j\hfill \end{array}$.

However, due to the non-linear unfolding rate function and forces we cannot solve Eqs (9) and (10) analytically and the numerical solution is unstable. Therefore, it is already highly costly to compute *P*_{1}, and costs are rising with every unfolded Ig domain. We therefore define basic probabilities which add up to the desired ones. Let *p*_{i1i2….im}(*t*) be the probability that at a given time *t*, *i*_{1} proteins of the first cluster, *i*_{2} of the second, …., and *i*_{m} of the *m*-th cluster are unfolded. The dynamics of these probabilities can be described with a simple system of linear ordinary differential Eq (11).

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill d{p}_{{i}_{1}....{i}_{m}}/dt& =& -\left(\sum _{l=1}^{m}({n}_{l}-{i}_{l})\xb7{u}_{l}\left({F}_{{i}_{1}+...+{i}_{m}}\left(t\right)\right)\right){p}_{{i}_{1}....{i}_{m}}(1-{\delta}_{{i}_{l},{n}_{l}})\hfill \\ & & +({n}_{1}-{i}_{1}+1){u}_{1}\left({F}_{{i}_{1}+...+{i}_{m}-1}\left(t\right)\right){p}_{({i}_{1}\phantom{\rule{-0.166667em}{0ex}}-\phantom{\rule{-0.166667em}{0ex}}1){i}_{2}...{i}_{m}}(1-{\delta}_{{i}_{1},0})+...\hfill \\ & & +({n}_{m}-{i}_{m}+1){u}_{1}\left({F}_{{i}_{1}+...+{i}_{m}-1}\left(t\right)\right){p}_{{i}_{1}{i}_{2}...({i}_{m}\phantom{\rule{-0.166667em}{0ex}}-\phantom{\rule{-0.166667em}{0ex}}1)}\left(t\right)(1-{\delta}_{{i}_{m},0}),\hfill \\ & & 0\le {i}_{1}\le {n}_{1},...,0\le {i}_{m}\le {n}_{m}\hfill \end{array}\end{array}$$

(11)

This detour has the advantage that we can solve the system fast and stable by an implicit Euler method. Finally, we calculate

$$\begin{array}{c}\hfill {P}_{l}=\sum _{\lambda \in {\Lambda}_{l}}{p}_{\lambda},\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}l=0,...,N,\end{array}$$

(12)

where ${\Lambda}_{l}=\{({i}_{1},...,{i}_{m}),0\le {i}_{1}\le {n}_{1},...,0\le {i}_{m}\le {n}_{m},{\sum}_{k=1}^{m}{i}_{k}=l\}$. The expected force in a half sarcomere is then

$$\begin{array}{c}\hfill E\left(F\left(t\right)\right)=\sum _{l=0}^{N}{P}_{l}\left(t\right){F}_{l}\left(t\right).\end{array}$$

(13)

Again, if refolding has to be taken into account the adjustment of the formulae is straight forward. For example, the formulation of refolding and folding of one cluster reads

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill d{p}_{i}/dt& =& (n-i+1)u\left({F}_{i-1}\left(t\right)\right){p}_{i-1}\left(t\right)+(n-i-1)r\left({F}_{i+1}\left(t\right)\right){p}_{i+1}\left(t\right)\hfill \\ & & -(n-i)u\left({F}_{i}\left(t\right)\right){p}_{i}\left(t\right)-ir\left({F}_{i}\left(t\right)\right){p}_{i}\left(t\right),\hfill \end{array}\end{array}$$

(14)

where *p*_{i}(*t*) is the probability that *i* Ig domains are folded at time *t* and *r* is the refolding rate.

Briefly, single myofibrils isolated from rabbit psoas muscle were wrapped around a glass needle at one end and fixed to a nanolever at the other end. The experimental set-up allowed for length change and force measurements, respectively. Myofibrils were passively stretched in a non-activating (pCa 8.0) solution containing ATP. The experimental procedure is described in detail elsewhere [25, 43]. Since extracellular connective tissues are absent in a single myofibril, all forces are attributed to proteins comprising sarcomeres in series. Due to sarcomere non-uniformities the scaling from a half sarcomere to a myofibril is not trivial. However, properties like hysteresis or the smoothing of unfolding events in a half sarcomere are reflected in single myofibrils [18]. Ethics approval for the single myofibril experiments was granted by the Life and Environmental Sciences Animal Ethics Committee of the University of Calgary.

The first simulations are based on 5 different unfolding clusters, each cluster covering 10 proximal Ig domains. Monte Carlo simulations of single titin strands show the typical saw-tooth pattern [5, 9, 36, 38] while this pattern is completely averaged out in a short myofibril comprised of 7 sarcomere (Fig 3).

The comparison between Monte Carlo approximation based on the simulation of 200 titin strands and the exact solution reveals that mean forces of the Monte Carlo simulation and the exact solution correspond almost perfectly. In fact, the difference between the mean force of 200 titin strands and the exact solution stays within a range of 2pN for stretches from 1 to 2μm half sarcomere length. The comparison of the unfolding probabilities as a function of half sarcomere length and the corresponding histograms (bin size based on the Freedman-Diaconis rule [44]) reveal that the shape of the exact solution is preserved by Monte Carlo simulations (Fig 4).

While the characteristics of the expected forces are comparable between Monte Carlo approximation and the exact solution other characteristics, like the most likely force at the first unfolding event, deviate (Fig 5) demonstrating that the simulation of more than 200 titin strands is necessary to capture the details of complex unfolding dynamics of all titin strands in a half sarcomere.

Furthermore, the complex folding and refolding behaviour of Ig domains is captured with both methods. As an example, we simulated the following experiment which was based on a study on single titin molecules reported in [5] (Fig 6): A half sarcomere is stretched from 1μm to 1.7μm half sarcomere length and re-shortened from 1.7μm to 1μm in two subsequent cycles. After a rest period of 30 seconds a third stretch-relaxation cycle was performed. Hysteresis declined in the second compared to the first cycle, but was fully recovered in the third cycle following a 30s rest; in single titin molecule experiments, the hysteresis in the third cycle even exceeded the hysteresis observed in the first cycle. Its worth pointing out that, unlike single titin molecule experiments [5], we cannot predict higher forces in any cycle following the first one with our exact solution, as we assume that all Ig domains are folded at the start of our simulations (Fig 6).

Due to the architecture of a single myofibril, where a (possibly large) number of half sarcomeres are aligned in series, the comparison between simulations on the half sarcomere level and experimental results on the myofibril level have to be done carefully. However, since the main characteristics of the simulations should be preserved on the myofibrillar level, the comparison gives us the opportunity to compare model predictions to the passive behaviour of sarcomeres in their natural environment. As an example, we stretched five single myofibrils to an average sarcomere length of approximately 4μm and then immediately performed 10 stretch-shortening cycles of an amplitude of around 0.25μm. We can observe a typical decline in peak forces over the 10 stretch-shortening cycles [45]. The corresponding Monte Carlo simulations (based on 200 titin strands) show a comparable qualitative decline of peak forces; the results differ quantitatively as we do not take sarcomere non-uniformities into account. Moreover, we presumed a simplified length behavior to prevent artefacts (see Fig 7).

We were interested in an efficient algorithm to compute average passive forces in half sarcomeres. We analysed the model by seeking to recover primary characteristics of experiments of single titin molecules as well as single myofibrils rather than comparing absolute values. We neither performed parameter estimation, nor did we try to fit experimental data. However, in terms of parameter estimation, the proposed algorithm should prove a highly useful alternative to classic Monte Carlo simulations.

There are certain limitations to the proposed algorithm for determining the exact solution. When the number of Ig clusters is high, Monte Carlo simulations might be the faster option. On the other hand, the smaller the number of clusters, the faster is the calculation of the exact solution compared to the simulation of hundreds of titin strands. In the extreme case, when each Ig domain is attributed a significantly different mechanical stability state and allowed to refold, the exact solution is numerical expensive for long stretches. Therefore, Monte Carlo simulations are the preferred approach. On the other hand, if one is interested in the probability distributions of unfolding forces of the first unfolding events, the exact solution still provides an efficient framework. The same is true for parameter estimations.

We introduce domain clusters in order to simplify the complexity of the folding / refolding process. While a limited number of clusters is favourable in terms of calculation time, it might be an over simplification. To study whether this simplification has a significant impact on the simulation outcome, we compare average forces and the probability of the first three unfolding events for two different cluster models. First, we look at 50 different clusters with linearly declining spontaneous unfolding rate constants and a random Gaussian perturbation and compare the results to simulations based on 5 different clusters with equidistant spontaneous unfolding rate constants. Since the overall shape of spontaneous unfolding rate constants are fairly similar for the two simulations, the corresponding average forces and unfolding forces are also remarkably similar (Fig 8).

However, clusters have to be chosen carefully. While the choice of the right clusters is of less importance when looking at qualitative behaviour of our system—we could perfectly explain hystereses of single titin molecules with one cluster—it is crucial in terms of parameter estimation or quantitative analysis. A perfect example would be the low force unfolding reported in a recent study by Martonfalvi et al [5]. The unfolding rates and first Ig clusters have to be chosen accordingly in order to cover the whole range of unfolding behavior that has been observed experimentally.

Finally, we would like to point out that the framework presented in this work is independent of the force model and the unfolding / folding rate functions used. Our approach is applicable whenever the average force of complex proteins is of interest.

Funding provided by the Austrian Science Fund (FWF), Grant Number T478-N13: http://www.fwf.ac.at, Canada Research Chair: http://www.chairs-chaires.gc.ca, Killam Foundation: http://killamlaureates.ca/, Natural Sciences and Engineering Research Council of Canada: http://www.nserc-crsng.gc.ca, and the Canadian Institutes of Health Research: http://www.cihr-irsc.gc.ca. The funders had no role in study design, data collection and analysis or decision to publish.

Data Availability

All relevant data are within the paper.

1.
Granzier H, Labeit S. Structure-function relations of the giant elastic protein titin in striated and smooth muscle cells. Muscle Nerve. 2007.
December;36(6):740–755. doi: 10.1002/mus.20886
[PubMed]

2.
Linke W. Stretching molecular springs: elasticity of titin filaments in vertebrate striated muscle. Histol Histopathol. 2000.
July;15(3):799–811. [PubMed]

3.
Bennett P, Hodkin T, Hawkins C. Evidence that the tandem Ig domains near the end of the muscle thick filament form an inelastic part of the I-band titin. J Struct Biol. 1997;120(1):93–104. doi: 10.1006/jsbi.1997.3898
[PubMed]

4.
Trinick J. End-filaments—A new structural element of vertebrate skeletal-muscle thick filaments. J Mol Biol. 1981;151(2):309–314. doi: 10.1016/0022-2836(81)90517-9
[PubMed]

5.
Martonfalvi Z, Bianco P, Linari M, Caremani M, Nagy A, Lombardi V, et al.
Low-force transitions in single titin molecules reflect a memory of contractile history. J Cell Sci. 2014.
February
15;127(4):858–870. doi: 10.1242/jcs.138461
[PubMed]

6.
Rief M, Gautel M, Oesterhelt F, Fernandez J, Gaub H. Reversible unfolding of individual titin immunoglobulin domains by AFM. Science. 1997.
May
16;276(5315):1109–1112. doi: 10.1126/science.276.5315.1109
[PubMed]

7.
Carrion-Vazquez M, Oberhauser A, Fowler S, Marszalek P, Broedel S, Clarke J, et al.
Mechanical and chemical unfolding of a single protein: A comparison. Proc Natl Acad Sci U S A. 1999.
March
30;96(7):3694–3699. doi: 10.1073/pnas.96.7.3694
[PubMed]

8.
Tskhovrebova L, Trinick J, Sleep J, Simmons R. Elasticity and unfolding of single molecules of the giant muscle protein titin. Nature. 1997.
May
15;387(6630):308–312. doi: 10.1038/387308a0
[PubMed]

9.
Kellermayer M, Smith S, Granzier H, Bustamante C. Folding-unfolding transitions in single titin molecules characterized with laser tweezers. Science. 1997.
May
16;276(5315):1112–1116. doi: 10.1126/science.276.5315.1112
[PubMed]

10.
Marszalek P, Lu H, Li H, Carrion-Vazquez M, Oberhauser A, Schulten K, et al.
Mechanical unfolding intermediates in titin modules. Nature. 1999.
November
4;402(6757):100–103. doi: 10.1038/47083
[PubMed]

11.
Levitt M, Warhsel A. Computer-simulation of protein folding. Nature. 1975;253(5494):694–698. doi: 10.1038/253694a0
[PubMed]

12.
Lu H, Isralewitz B, Krammer A, Vogel V, Schulten K. Unfolding of titin immunoglobulin domains by steered molecular dynamics simulation. Biophys J. 1998.
August;75(2):662–671. doi: 10.1016/S0006-3495(98)77556-3
[PubMed]

13.
Best R, Fowler S, Herrera J, Steward A, Paci E, Clarke J. Mechanical unfolding of a titin Ig domain: Structure of transition state revealed by combining atomic force microscopy, protein engineering and molecular dynamics simulations. J Mol Biol. 2003.
July;330(4):867–877. doi: 10.1016/S0022-2836(03)00618-1
[PubMed]

14.
Gao M, Lu H, Schulten K. Unfolding of titin domains studied by molecular dynamics simulations. J Muscle Res Cell Motil. 2002;23(5–6):513–521. doi: 10.1023/A:1023466608163
[PubMed]

15.
Higuchi H, Nakauchi Y, Maruyama K, Fujime S. Characterization of beta-connectin (titin-2) from striated-muscle by dynamic ligth-scattering. Biophys J. 1993.
November;65(5):1906–1915. doi: 10.1016/S0006-3495(93)81261-X
[PubMed]

16.
Linke W, Ivemeyer M, Mundel P, Stockmeier M, Kolmerer B. Nature of PEVK-titin elasticity in skeletal muscle. Proc Natl Acad Sci U S A. 1998.
July
7;95(14):8052–8057. doi: 10.1073/pnas.95.14.8052
[PubMed]

17.
Minajeva A, Kulke M, Fernandez J, Linke W. Unfolding of titin domains explains the viscoelastic behavior of skeletal myofibrils. Biophys J. 2001.
March;80(3):1442–1451. doi: 10.1016/S0006-3495(01)76116-4
[PubMed]

18.
Herzog JA, Leonard TR, Jinha A, Herzog W. Are titin properties reflected in single myofibrils?
J Biomech. 2012.
July
26;45(11):1893–1899. doi: 10.1016/j.jbiomech.2012.05.021
[PubMed]

19.
Fernandez J, Li H. Force-clamp spectroscopy monitors the folding trajectory of a single protein. Science. 2004.
March
12;303(5664):1674–1678. doi: 10.1126/science.1092497
[PubMed]

20.
Oberhauser A, Hansma P, Carrion-Vazquez M, Fernandez J. Stepwise unfolding of titin under force-clamp atomic force microscopy. Proc Natl Acad Sci U S A. 2001.
January
16;98(2):468–472. doi: 10.1073/pnas.98.2.468
[PubMed]

21.
Popa I, Kosuri P, Alegre-Cebollada J, Garcia-Manyes S, Fernandez JM. Force dependency of biochemical reactions measured by single-molecule force-clamp spectroscopy. Nat Protoc. 2013.
July;8(7):1261–1276. doi: 10.1038/nprot.2013.056
[PMC free article] [PubMed]

22.
Brujic J, Hermans RIZ, Garcia-Manyes S, Walther KA, Fernandez JM. Dwell-time distribution analysis of polyprotein unfolding using force-clamp spectroscopy. Biophys J. 2007.
April;92(8):2896–2903. doi: 10.1529/biophysj.106.099481
[PubMed]

23.
Zhmurov A, Dima RI, Barsegov V. Order Statistics Theory of Unfolding of Multimeric Proteins. Biophys J. 2010.
September
22;99(6):1959–1968. doi: 10.1016/j.bpj.2010.07.012
[PubMed]

24.
Joumaa V, Leonard TR, Herzog W. Residual force enhancement in myofibrils and sarcomeres. Proc R Soc Lond B Biol Sci. 2008.
June
22;275(1641):1411–1419. doi: 10.1098/rspb.2008.0142 [PMC free article] [PubMed]

25.
Leonard TR, Herzog W. Regulation of muscle force in the absence of actin-myosin-based cross-bridge interaction. Am J Physiol Cell Physiol. 2010.
July;299(1):C14–C20. doi: 10.1152/ajpcell.00049.2010
[PubMed]

26.
Rassier D, Herzog W, Pollack G. Dynamics of individual sarcomeres during and after stretch in activated single myofibrils. Proc R Soc Lond B Biol Sci. 2003.
August
22;270(1525):1735–1740. doi: 10.1098/rspb.2003.2418 [PMC free article] [PubMed]

27.
Minozzo FC, Baroni BM, Correa JA, Vaz MA, Rassier DE. Force produced after stretch in sarcomeres and half-sarcomeres isolated from skeletal muscles. Sci Rep. 2013.
July
31;3
doi: 10.1038/srep02320
[PMC free article] [PubMed]

28.
Rassier DE. Pre-power stroke cross bridges contribute to force during stretch of skeletal muscle myofibrils. Proc R Soc Lond B Biol Sci. 2008.
November
22;275(1651):2577–2586. doi: 10.1098/rspb.2008.0719 [PMC free article] [PubMed]

29.
Schappacher-Tilp G, Leonard T, Desch G, Herzog W. A Novel Three-Filament Model of Force Generation in Eccentric Contraction of Skeletal Muscles. PLoS One. 2015.
March
27;10(3). doi: 10.1371/journal.pone.0117634
[PMC free article] [PubMed]

30.
Campbell SG, Hatfield PC, Campbell KS. A Mathematical Model of Muscle Containing Heterogeneous Half-Sarcomeres Exhibits Residual Force Enhancement. PLoS Comput Biol. 2011.
September;7(9). doi: 10.1371/journal.pcbi.1002156
[PMC free article] [PubMed]

31.
Campbell KS. Interactions between Connected Half-Sarcomeres Produce Emergent Mechanical Behavior in a Mathematical Model of Muscle. PLoS Comput Biol. 2009.
November;5(11). doi: 10.1371/journal.pcbi.1000560
[PMC free article] [PubMed]

32.
Nishikawa KC, Monroy JA, Uyeno TE, Yeo SH, Pai DK, Lindstedt SL. Is titin a ‘winding filament’? A new twist on muscle contraction. Proc R Soc Lond B Biol Sci. 2012.
March
7;279(1730):981–990. doi: 10.1098/rspb.2011.1304 [PMC free article] [PubMed]

33.
Herzog W. The role of titin in eccentric muscle contraction [Editorial Material]. J Exp Biol. 2014.
August;217(16):2825–2833. doi: 10.1242/jeb.099127
[PubMed]

34.
Freiburg A, Trombitas K, Hell W, Cazorla O, Fougerousse F, Centner T, et al.
Series of exon-skipping events in the elastic spring region of titin as the structural basis for myofibrillar elastic diversity. Circ Res. 2000.
June
9;86(11):1114–1121. doi: 10.1161/01.RES.86.11.1114
[PubMed]

35.
Rief M, Fernandez J, Gaub H. Elastically coupled two-level systems as a model for biopolymer extensibility. Phys Rev Lett. 1998.
November
23;81(21):4764–4767. doi: 10.1103/PhysRevLett.81.4764

36.
Li H, Linke W, Oberhauser A, Carrion-Vazquez M, Kerkviliet J, Lu H, et al.
Reverse engineering of the giant muscle protein titin. Nature. 2002.
August
29;418(6901):998–1002. doi: 10.1038/nature00938
[PubMed]

37.
Makarov D, Hansma P, Metiu H. Kinetic Monte Carlo simulation of titin unfolding. J Chem Phys. 2001.
June
1;114(21):9663–9673. doi: 10.1063/1.1369622

38.
Bianco P, Martonfalvi Z, Naftz K, Koszegi D, Kellermayer M. Titin Domains Progressively Unfolded by Force Are Homogenously Distributed along the Molecule. Biophys J. 2015.
July
21;109(2):340–345. doi: 10.1016/j.bpj.2015.06.002
[PubMed]

39.
Marko J, Siggia E. Stretching DNA. Macromolecules. 1995.
December
18;28(26):8759–8770. doi: 10.1021/ma00130a008

40.
Bustamante C, Marko J, Siggia E, Smith S. Entropic elasticity of lambda-phage DNA. Science. 1994.
September
9;265(5178):1599–1600. doi: 10.1126/science.8079175
[PubMed]

41.
Leake M, Wilson D, Gautel M, Simmons R. The elasticity of single titin molecules using a two-bead optical tweezers assay. Biophys J. 2004.
August;87(2):1112–1135. doi: 10.1529/biophysj.103.033571
[PubMed]

42.
Daniel T, Trimble A, Chase P. Compliant realignment of binding sites in muscle: Transient behavior and mechanical tuning. Biophys J. 1998.
April;74(4):1611–1621. doi: 10.1016/S0006-3495(98)77875-0
[PubMed]

43.
Powers K, Schappacher-Tilp G, Jinha A, Leonard T, Nishikawa K, Herzog W. Titin fore is enhanced in actively stretched skeletal muscle. JEB. 2014;217:3629–3636. doi: 10.1242/jeb.105361 [PubMed]

44.
Freedman D, Diaconis P. On the histogram as a density estimator—L2 theory. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete. 1981;57(4):453–476. doi: 10.1007/BF01025868

45. Herzog J, Leonard TR, Jinha A, Herzog W; Biophys Soc. Titin Hysteresis is Greater for Actively Lengthened Compared to Passively Lengthened Skeletal Muscle Sarcomeres [Meeting Abstract]. Biophys J. 2015 Jan 27;108(2, 1):460A. 59th Annual Meeting of the Biophysical-Society, Baltimore, MD, Feb 07–11, 2015.

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |