Resolution and uniqueness of estimated parameters of a model of thin filament regulation in solution
The estimation of chemical kinetic rate constants for any non-trivial model is complex due to the nonlinear effects of second order chemical reactions. We developed an algorithm to accomplish this goal based on the Damped Least Squares (DLS) inversion method and then tested the effectiveness of this method on the McKillop-Geeves (MG) model of thin filament regulation. The kinetics of MG model is defined by a set of nonlinear ordinary differential equations (ODEs) that predict the evolution of troponin-tropomyosin-actin and actin-myosin states. The values of the rate constants are estimated by integrating these ODEs numerically and fitting them to a series of stopped-flow pyrene fluorescence transients of myosin-S1 fragment binding to regulated actin in solution. The accuracy and robustness of the estimated rate constants are evaluated for DLS and two other methods, namely quasi-Newton (QN) and simulated annealing (SA). The comparison of these methods revealed that SA provides the best estimates of the model parameters because of its global optimization scheme. However it converges slowly and does quantify the uniqueness of the estimated parameters. On the other hand the QN method converges rapidly but only if the initial guess of the parameters is close to the optimum values, otherwise it diverges. Overall, the DLS method proves to be the most convenient method. It converges fast and was able to provide excellent estimates of kinetic parameters. Furthermore, DLS provides the model resolution matrix, which quantifies the interdependence of model parameters thereby evaluating the uniqueness of their estimated values. This property is essential for estimating of the dependence of the model parameters on experimental conditions (e.g. Ca2+ concentration) when it is assessed from noisy experimental data such as pyrene fluorescence from stopped-flow transients. The advantages of the DLS method observed in this study should be further examined in other physicochemical systems to firmly establish the observed effectiveness of DSL vs. the other parameter estimation methods.
1. Introduction<!>2. Methods<!>2.1. The McKillop-Geeves three-state model of thin filament regulation<!>2.2.1. Damped Least Squares (DLS)<!>2.2.1.1. The resolution matrix (RM)<!>2.2.2. Quasi-Newton (QN)<!>2.2.3. Simulated annealing (SA)<!>2.2.4. Convergence criteria<!>2.2.5. Effective kinetic binding rate constant, kobs<!>3. Results<!>3.1. Estimation of MG model parameters by DLS, quasi-Newton and SA estimation methods<!>3.2. Calcium concentration dependence of KB and KT<!>3.3. Estimation errors and convergence of the estimated parameters<!>4. Discussion and conclusions
<p>The kinetics of a chemically reacting system is usually modelled using ordinary differential equations that are parameterized by a reduced set of reaction rate constants. Precise knowledge of these constants is required to characterize the dynamical behavior of the system. The estimation of the reaction rate constants from experimental data constitutes an inverse problem and is much more complicated than solving the model for a given set of rate constants. Inverse problems arising in chemical kinetics can be addressed by discrete inverse theory (Menke, 1989). The most common methods for parameter estimation on chemical kinetics are least squares methods that minimize an objective function iteratively by working with its gradient (Farinha et al., 1997; Lisy and Simon, 1998; Tadi and Yetter, 1998). In these methods the objective function is defined as the error between the experimental observations and the predictions of the model, which depends on the model parameters. Several other approaches have been developed to solve this difficult minimization problem. One of them is the so-called "trust region", which uses successive quadratic programming (SQP), as proposed by Arora and Biegler (2004). Also in use are various stochastic search techniques, including simulated annealing (Eftaxias et al., 2002; Goffe et al., 1994) and genetic algorithms (Terry and Messina, 1998).</p><p>Ideally, the parameter set that minimizes the objective function is a manifold of dimension zero in which the minimum of error is a single distinguished point in the error landscape. In this case, the error increases as the parameter values vector moves away from this point in any arbitrary direction. However, in systems where multiple parameters need to be estimated, the optimal parameter set may span a manifold of dimensions equal or greater than one, forming "valleys" or "hyper-valleys" in the error landscape. As a result, the model output may be independent of changes in one parameter or a certain combination of parameters, which prevents the univocal estimation of these parameters. A sensitivity analysis is therefore necessary to explore the variations in the model output with perturbations in the parameters, which are defined by the sensitivity coefficients (Frank, 1978; Tomovic and Vukobratovic, 1972; Varma et al., 1999). In fact, the parameters of the model can be determined only if the sensitivity coefficients are non-zero and linearly independent (Beck and Arnold, 1977). More rigorously, Tang et al. (2005) used the "physically bounded Gauss-Newton" (PGN) method by utilizing the sensitivity matrix as a global map to determine the unknown kinetic parameters of complex reaction networks that contain dozens of species and hundreds of reactions.</p><p>The primary focus of this study is on the iterative Damped Least Squares (DLS) method which estimates kinetic parameters and provides a model resolution matrix (RM) (Menke, 1989). The DLS method, which is also known as the Levenberg-Marquardt method (Kecman, 2001; Nocedal and Wright, 2006), utilizes the sensitivity matrix as a search map in the same way as the PGN method, and provides information about the parameter uniqueness via the model RM. We applied the DLS method to determine the parameters of a kinetic model by fitting experimental data from the isotherm of the myosin-S1 binding to regulated F-actin in solution. The fidelity of the DLS method was evaluated by comparing its results with those obtained from two other algorithms: simulated annealing, and quasi-Newton. The predictions of the kinetics model are calculated by a probabilistic algorithm based on the McKillop-Geeves (MG) three-state model of thin filament regulation in solution (McKillop and Geeves, 1993). According to this model, the actin-associated regulatory protein complexes consisting of tropomyosin (Tm) and troponin (Tn) switch between the three azimuthal positions on the actin filament surface thereby preventing or allowing myosin-S1 to (weakly or strongly) bind to actin. The kinetics of the interconversion of TmTn-actin states depends on calcium concentration and represents the primary mechanism of regulation of contraction in striated muscles. The estimated key parameters of the MG model and their dependence on calcium by the three-parameter methods are compared, and their uniqueness and robustness is discussed. The rate of convergence of the estimated parameters and the effectiveness of different parameter estimation methods is also evaluated.</p><!><p>All experimental data contain experimental error and natural fluctuations. In order to better understand how and to what degree these uncertainties affect estimated model parameters, we applied a probabilistic algorithm for the simulation of the MG three-state model of thin filament regulation in solution (McKillop and Geeves, 1993) developed by Chen et al. (2001). We tested the robustness of DLS parameter estimation by comparing the model predictions with two other widely used parameter estimation methods, namely, the quasi-Newton and simulated annealing. The model predictions required by the three estimation methods were calculated using the probabilistic formulation of the MG model (Chen et al., 2001). This method was preferred because of its speed since many model calculations may be needed to reach convergence.</p><!><p>In vertebrate skeletal and cardiac muscles the interaction between myosin and actin is regulated by the actin-associated proteins, tropomyosin (Tm) and troponin (Tn), depending on the concentration of calcium (Ca2+). The soluble fragment of the myosin molecule S1 is widely used for studying kinetics of myosin binding to regulated actin filaments in solution. This fragment, also known as the motor domain, contains all of the ATPase and actin binding properties of the parent myosin. In the absence of nucleotide, myosin-S1 forms a tight (rigor-like) bond to actin filaments. McKillop and Geeves (1993) proposed that the regulation of tropomyosin and troponin-containing thin filaments can be interpreted using three-states of the actin filament: (1) the "blocked" state, in which myosin-S1 binding to actin is prohibited; (2) the "closed" state, in which S1 can bind with actin, but cannot be isomerized further to next step; (3) and in the "open" state, where no limitation to S1 binding to actin is imposed. In this study, the unit size of Tm-Tn complex is assumed to cover 7 actin monomers,1 denoted as actin7·TmTn (Maytum et al., 1999) (see Fig. 1). The repeat of TmTn every 7 units uniquely defines a TmTn unit that can rigidly move between the three states. Because Ca2+ binding to Tn significantly decreases the affinity of Tn to actin, the distribution between these three states is therefore affected by Ca2+ concentration. There are three myosin states in the model: one unbound state where myosin is in solution and two bound states where myosin is bound to actin. The bound states are denoted as weakly bound A-states and strongly bound, i.e. rigor-like states, as R-states (Fig. 1).</p><p>The MG model (Fig. 1) can be fully described by defining each state as combination of either, blocked, closed or open TmTn state, and the number actin sites with myosin bound in A- and R-stated within an actin7·TmTn unit. The complete chemical kinetics of the MG model can be described by 45 states (see Fig. 2 for details) and 45 corresponding chemical kinetics equations, as previously described by Chen et al. (2001). These equations are
(1)dp1(t)dt=−k+Bp1+k−Bp2,dp2(t)dt=k+Bp1−(k−B+7k+1c+k+T)p2+k−1p3+k−Tp10,⋮dp45(t)dt=k+2p44−7k−2p45, where pi(t) is the fraction TmTn units in state i, which is defined by position of Tm (i.e. in block, closed or open state) and a particular combination of actin unoccupied sites and S1 bound in A- or R-states within the actin7·TmTn unit. The equilibrium constants of the model are defined as KB = k+B/k−B, KT = k+T/k−T, K1 = k+1/k−1, and K2 = k+2/k−2, where k+B, k+T, k+1, and k+2 are forward rate constants, and k−B, k−T, k−1, and k−2 are backward rate constants. The system of Eq. (1) is nonlinear because the forward transition rate between the free actin state and myosin bound to actin in the Astate depends on the concentration of unbound S1. In Eq. (1) we denoted S1 concentration as c = [S1], thus the effective myosin binding rate constant, k+1c, is in units s−1. The concentration of S1 in solution, i.e. the concentration of unbound myosin, decreases as myosin binds to actin and, therefore, decreases the effective rate of myosin binding. The concentration of free myosin-S1 in solution is equal to the initial concentration of free myosin-S1, [S1O], minus concentrations of myosin to bound to actin in either A- or in R-state. Thus the binding rate depends on probabilities of all myosin bound states, pbound. The equations of the MG model are solved numerically by Gear's backward differentiation formulas (up to order five) (Hindmarsh, 1972). All 45 equations can be easily solved for a wide range of parameter combinations. However, in some extreme cases in which some of the parameters take much larger values than others, the set of Eq. (1) becomes numerically "stiff" and it is difficult to solve using standard methods. In those cases, Monte Carlo simulations should be used in order to solve the resulting stiff set of equations efficiently.</p><p>Once we determine the vector of 45 states, p(k, t), where vector k = (k1,…kn) represents an array of rate transition constants and other relevant model parameters, we calculate the fraction of actin sites in each one of the three actomyosin states by summing the occupancy of the actin sites over all TmTn units, and normalizing by the total number of actin sites occupied by S1 in the R-state when the system is in equilibrium. The model predictions are tested against measurements of the pyrene fluorescence intensity during stopped-flow experiments (Boussouf et al., 2007a,b; Boussouf and Geeves, 2007). A drop in pyrene fluorescence is proportional to myosin binding to actin in the R-state. Thus, the calculated instantaneous fractions of actin sites that are not occupied or which are in the A-state (i.e. which are not in the R-state), denoted g(κ, t), can be compared with corresponding experimental data, dobs(t), at the same instant. Here, the vector κ = (κ1, … κm) represents the set of m free model parameters that need to be estimated. Note that κ is usually a subset of the complete parameter set k, since some of the parameters are prescribed by the experimental protocol, in which case m ≤ n. For example, the prescribed concentration of actin, myosin and calcium are the same in both the model simulations and those used in the experiments. Also, m is reduced if some parameters, such as several rate or equilibrium state transition constants, vary a little over the course of multiple experiments and can be measured independently.</p><!><p>Damped Least Square (i.e. Levenberg–Marquardt) inversion is a widely used method for the estimation of model parameters that has two important features: quantitative evaluation of the uniqueness of the estimated parameters and good parameter resolution (Menke, 1989). The DLS method is based on the iterative minimization of the mean-square error of the model predictions with respect to experimental observations. In this study, we estimate the rate transition constants of the McKillop-Geeves model (Fig. 1) by minimizing the variance between predicted history of the fraction of actin sites unoccupied by bound myosin in R-state, g(κ, t), and the same fraction deduced from a time course of (normalized) pyrene fluorescence intensity, dobs(t). We exclusively fitted data recorded during stopped-flow experiments for various concentrations of Ca2+ in which actin concentration is in excess of myosin-S1 concentration (Boussouf et al., 2007a,b; Boussouf and Geeves, 2007). The measurements performed in these experiments covered a range of independent variables wide enough to uniquely resolve the set of model parameters for each Ca2+ concentration. These observations dobs(t) are represented by a set on N values of the pyrene fluorescence sampled at discrete time points.</p><p>The model predictions deviate from the experimental observations by e(κ, t) = dobs(t) − g(κ, t), yielding the mean-square-error
(2)E(κ)=‖e‖22N=1N∑i=1N(gi(κ)−diobs)2, which represents an integral measure of the "goodness" of the model prediction fit to the observations. Minimization of Eq. (2) provides the set of free parameters κ that best fit the experimental data. Finding the best set of parameters which minimize Eq. (2) by the DLS method leads to the following iterative scheme
(3)Δκr+1=[GrTGr+ε2I]−1GrTer=Jerκr+1=κr−Δκr+1 where er is the model error for the set of parameters obtained in the rth iteration, κr. Gr is the Jacobian matrix of the model predictions with respect to the parameters at iteration r,(Gr)i,j = ∂g(κr, ti)/∂κj; δκr+1 is the vector of estimated increments of parameter κr, and κr+1 the vector of estimated parameter values at iteration r +1. At the end of the iteration the error is calculated for the set of current parameters, κr, which is represented in the error landscape as a step towards the location of the minimal error. The r + 1th step is determined by the product of the error vector at the rth step, er and the pseudoinverse matrix
(4)Jr=[GrTGr+ε2I]−1GrT=Mr−1GrT, which depends on the local topology of the error landscape and provides a search "map" or solution direction in this landscape (see Appendix A for more details). The parameter ε2 is called the damping parameter and can take an arbitrary small value. This parameter ensures the invertability of Mr by keeping its smallest eigenvalue equal to, or greater than, ε2. In other words, ε2 ensures the stability of the iterative procedure by limiting the length of each iteration step in the error landscape to 1/ε2. As a counterpart, non-zero values of ε may slow convergence and/or degrade the resolution of the parameters because the objective function that is actually minimized for ε >0 is Φ = E(κ) + ε2|| Δκ||, instead of just E(κ). A suitable small value of ε should "damp" any disturbance (i.e. experimental data noise) and only slightly affect the model parameter values.</p><p>When the parameter values, and therefore the increments Δκ, vary in a wide range it is more efficient to use logarithmic variations (Levenberg, 1944). In this approach we look for the solution of log(g(κ, t)) ≅ log(dobs(t)) and the model deviation from the observations is e(κ, t) = log(dobs(t)/g(κ, t)). Also, the partial derivatives that define Gr are of log(g(κ, t)) with respect to log(κ), and the new value of the estimated parameter is κr+1 = κr exp(−Δ κr+1/κr) (see Appendix A for details).</p><!><p>The sensitivity of the model to each of the estimated parameters is examined by the resolution matrix R (an m × m matrix) at convergence: 
(5)Rr=JrGr=[GrTGr+ε2I]−1GrTGr</p><p>When the matrix 
GrTGr is a regular matrix and ε2 is small, then R ~ I and the estimated model parameters are uniquely determined, i.e. the obtained model is close to being correct. However, if the matrix R has appreciable non-zero off-diagonal elements, variations of a particular parameter can be compensated by adjustments of other parameters that fit the experimental data equally well. Furthermore, if all the elements in a row of the RM are small (|RM| ≪ 1), then the corresponding parameter cannot be estimated for the given data set because the error is affected very little, even by large changes in this parameter. In these two situations, the optimal set of parameters occupies a manifold of dimension ≥ 1 and forms a "valley" or "hypervalley" in the error landscape. The RM discriminates meaningless estimates and also indicates those parameters that may vary in a wide range without alternation of the system response, as one can move along a valley in the error landscape without modifying the mean-square error.</p><!><p>This method estimates the parameter set κ = (κ1, … κm) by minimizing the error E(κ, t) with the Newton algorithm (Dennis and Schnabel, 1983; Nocedal and Wright, 2006). In contrast with the standard Newton method, which requires numerical calculation of the gradient and Hessian of E(κ, t) in κ space, the QN method only requires computing the gradient; the Hessian is estimated using the gradients from successive iterations (Press et al., 1986). We used the implementation of the QN method provided by the Microsoft Power Fortran routine, BCONF, which includes an active set strategy (Gill and Murray, 1976). For a given error function, the QN method iteratively computes the search direction by using the step length determined by a local optimization of the function, called a line search, until it reaches the minimum of E(κ, t) function values. The Hessian is well-conditioned if the approximation by a quadratic function of the objective function surface E(κ, t) is the correct one which will be the case in the neighborhood of a local or global minimum when the convergence rate is quadratic (Dennis and Schnabel, 1989; Nocedal and Wright, 2006). Quite often, however, the objective function is poorly approximated by a quadratic function in the vicinity of the initial guess, especially if the latter is far from the optimum. Consequently, despite fast convergence, many trial and error runs involving different initial guesses may be necessary to determine the set of model parameters precisely.</p><!><p>This algorithm provides a global search method specifically designed for functions with multiple minima over wide range of parameters2 (Kirkpatrick et al., 1983; Press et al., 1986). This feature is especially important in chemical kinetics systems with many reaction equations and reaction rate constants which may vary considerably. Using a global strategy, each step of the SA algorithm replaces the current solution by a randomly chosen solution that is "nearby" in the parameter space; if the new solution is better it is chosen the new solution; whereas if it is worse the new solution can still be chosen, with a probability that depends on the difference between the corresponding function values and on a global parameter T (called the temperature) that is gradually decreased during the iterative process. As a result, the iterating solution changes almost randomly when T is large, but increasingly "down the well" as T goes to zero. The allowance for a random "uphill" enables escaping from local minima. A similar, although less powerful global search can be performed in DLS by gradually increasing the damping parameter ε2, as the iteration process proceeds. We used the SA implementation provided in the Microsoft Power Fortran routine, SA (Goffe et al., 1994) to search a global minimum of the error function E(κ, t). This powerful method is effective in finding an absolute minimum of E(κ, t) but requires a large number of model evaluations. Thus its application is limited to parameter estimation in models whose equations can be solved quickly.</p><!><p>In DLS, the iterative process is stopped when the error and the norm of the increment of the estimated parameters are less than the prescribed respective tolerances. Also the maximal number of iterations is usually limited (here, the limit was 100 and the convergence is typically reached in less than 20 steps). The proprietary QN and SA implementations have their own exit strategies. There are two stopping criteria for BCONF that occur, either when the norm of the gradient is less than a given gradient tolerance (<10−8), or when the scaled distance between the latest two steps is less than the step tolerance (<10−10). The program also stops when it reaches prescribed values for the maximum number of function evaluations (default: 100), or of gradient evaluations (default: 400). SA terminates when either one of these two exit conditions is satisfied: (1) when the error is less than default tolerance value of 10−6 or (2) when maximum number of function evaluations is reached (currently set to 105).</p><!><p>One important determinant of myosin-S1 binding to regulated thin filament via the TmTn complex is the effective kinetic binding rate constant, kobs. When the actin concentration is much larger than the myosin-S1 concentration ([A] ≫ [S1]) in the presence of calcium, a single exponential is observed and kobs = (−d[A]/[Ao])k+1[S1] − k−1[A.M] represents the rate of decrease of unbound myosin concentration (McKillop and Geeves, 1993). Here [A.M] is the concentration of actin sites in weakly bound S1, i.e., A-states. When KT is large (~200), the population of the A-state is small, and the reverse rate is negligible in absence of nucleotide, and thus the forward rate is approximately proportional to the rate of fluorescence decay in stopped-flow experiment (after subtraction of background fluorescence). Also the observed rate constant, kobs, is independent of the ratio between actin to S1 over the range 5:1 to 20:1. In this case for lower and intermediate concentrations of Ca2+, kobs can be calculated as</p><p>For the family of stopped-flow transient data, we can obtain the fractional 
kobsf by fitting the stopped-flow data. The fluorescence data are corrected for the background fluorescence then they are typically normalized by setting the equilibrium value to zero, followed by normalizing to the corrected initial fluorescence. Also for the effective comparison of 
kobsf from stopped-flow data with estimated 
kobsf from the best fits (6), kobs is further normalizing at any Ca2+ concentration to the highest Ca2+ concentration (i.e. pCa = 4.6) after subtracting the value of kobs at lowest concentration (i.e. pCa = 4.6): 
(7)kobsf=kobs(pCa)−kobs(8.9)kobs(4.6)−kobs(8.9).</p><p>This procedure allows a comparison of the fractional equivalent rate constant, 
kobsf, from the estimated MG model parameters by different estimation methods for every set of experimental kinetic data at a prescribed Ca2+ concentration.</p><!><p>To evaluate the effectiveness and robustness of the parameter estimation methods we fitted the model predictions to stopped-flow transients for S1 binding to pyrene-actin in the presence of regulatory proteins for a specified Ca2+ concentration (Boussouf et al., 2007a,b; Boussouf and Geeves, 2007). To improve accessibility of the parameter estimation methods for the MG model, we developed Visual Basic GUI, which integrates all computing components into a single shell application.3</p><!><p>The McKillop–Geeves (MG) three-state rigid chain model has four free parameters: the equilibrium constants KB, KT, K1, and K2 (see Section 2). We estimated these parameters using the DLS, QN and SA methods by fitting the predictions of the MG model to stopped-flow transients for binding S1 to excess actin ([A] = 2.5 μM, [S1] = 0.25 μM) (Boussouf et al., 2007a). The fitted curves for these three methods in the case of high or low Ca2+ differ little and agree well with the experimental data (Fig. 3). Table 1 shows the values of the four equilibrium constants obtained by DLS, QN and SA for high or low Ca2+ concentrations. It is important to note that different estimation methods estimated somewhat different of values of the parameters KB, KT, K1 and K2, although the differences in the fitting error are hardly noticeable (see Table 1) and the fits are almost identical (Fig. 3). The same behavior is observed at other Ca2+ concentrations too, thus a consistent trend of the parameter values with Ca2+ concentration is extremely difficult to establish if all four parameters are fitted simultaneously. The parameter values estimated from the different Ca2+ concentration by the above three estimation methods showed substantial differences in some cases, or did not show a consistent dependence on Ca2+ concentration.</p><p>This inconsistency may be attributed to a relatively flat error landscape, which would cause different sets of parameters to yield similar model predictions. When the differences between model predictions are within the experimental error (i.e. fluctuations in measured fluorescence), the parameter estimation methods would not able to discern which solution is better. In the SA paradigm, this effect can be interpreted as a residual temperature increment due to the entropy of the experimental measurements.</p><p>We explored the possibility of reducing this effect by filtering the experimental data through averaging pairs of two neighboring points. The difference between the parameters estimated from the original data and the filtered data was minor, the error was slightly reduced and the resolution matrix was improved in some cases. However, this apparent improvement was inconsistent, even when data were filtered a second time. Thus a reduction of experimental fluctuations has not resulted in significantly better convergence and better parameter resolution, pointing out that additional effects affected the goodness of the fit.</p><p>These findings indicate that the uniqueness of the sets of parameters estimated from stopped-flow transients should be evaluated carefully. The DLS method enables such an evaluation through analysis of the resolution matrices (see Table 2). Reducing the number of estimated parameters by, for example, fixing some model parameter values by using reliable values from other experiments (e.g. K1 and K2) improves the consistency of the remaining parameters with Ca2+ concentration. Fig. 4A shows that the best fits of stopped-flow transients obtained by the DLS, QN and SA methods at pCa of 8.9 are almost identical even when the number of estimated parameters is reduced to two. However, these fits are not as good as fits shown in Fig. 3 where all four free parameters are estimated. The reduction of the number of parameters from four to three and then to two also causes some minor differences in fits (Fig. 4B), but these differences are small and similar to the differences observed in Fig. 4A. Despite these differences, the reduction of the number of parameters brings the values of the estimated parameters closer by the different methods (Table 1), and improves the resolution matrix (Table 2), i.e. improves the robustness of the parameter estimation methods, but at the expense of less accurate fits (Fig. 4B, Table 1).</p><!><p>Establishing the dependence of the equilibrium constants KB and KT on the concentration of Ca2+ by estimating the parameters of the MG model from the fits to families of stopped-flow transient experiments is important for the quantitative understanding of thin filament regulation in intact muscle fibers. As discussed above, fixing K1 and K2 and estimating only KB and KT can improve the uniqueness of the estimated parameters. The rationale for this simplification is based on three observations made in Table 2 and Fig. 5: (i) K1 is uniquely determined even when fitting all four parameters of the MG model; (ii) K1 and K2 vary little with the Ca2+ concentration; and (iii) the best estimates from the best fits of the full set of parameters over eleven calcium concentrations are 0.25 μM−1 and 170, respectively. These values vary little with concentration and they are consistent with previous studies. McKillop and Geeves (1993) showed that the binding equilibrium constant, K1, from the kinetic data in the presence of Ca2+ at actin concentrations is between 0.1 and 0.2 μM−1. The same authors reported that the value of K2 is around 200 in the absence of nucleotide (McKillop and Geeves, 1993).</p><p>In order to illustrate how a reduction of the number of estimated parameters affects the values of the parameters and the quality of the fits we compared the values of the estimates of all four free parameters, of three free parameters with prescribed K1 = 0.25 μM−1, and of two free parameters with prescribed K2 = 170 (Table 1,K
Fig. 4B). Fixing 1 and K2, i.e. reducing the number of the free parameters to two is advantageous for global methods such as SA because it significantly reduces the number of model calculations and it enhances the uniqueness of the estimated KB, and KT over a wider range of Ca2+ concentrations (see DLS resolution matrices in Table 2 and Fig. 5). All parameters are estimated from stopped-flow data at actin concentrations in excess of S1 ([A] = 2.5 μM, [S1] = 0.25 μM).</p><p>The Kb values estimated by DLS, QN and SA display a sigmoidal Ca2+ dependence with the steepest variation at pCa between 6.6 and 5.4 (Fig. 6A). Hill's coefficients obtained from parameters obtained by each of three-parameter estimation methods are similar, averaging of about 2.25 (namely 2.15, 2.28 and 2.34 for DSL, SA and QN, respectively) and the 50% of KB range is achieved at pCa of 5.57, 5.58 and 5.64, respectively. The KT-Ca2+ dependence is less prominent, showing a mild increase from 0.03 at low Ca2+ concentration to ~0.125 at high Ca2+ (Fig. 6B). The relatively smooth trend is interrupted by two "spikes" at pCa = 6.0 and pCa = 5.4, which are likely caused by imperfections of the stopped-flow transients at early times, and the imprecise acquisition of the experiment starting time. These large variations of KT between neighboring concentrations and the variability of the KT values estimated by different methods are also reflected in the Hill coefficients, which are respectively, 0.47, 1.22 and 0.76 for DSL, SA and quasi-Newton. Half of the change in KT is achieved at lower Ca2+ concentrations (at pCa of 6.54, 6.54 and 6.27, respectively) than estimated from the KB-pCa relationship. These spikes and variations in the Hill coefficients, however, have almost no effect on the effective kinetic binding rate constant, kobs (Fig. 7).</p><p>The dependence of kobs with Ca2+ concentration is similar to the rate constant of single exponential decay of fractional fluorescence, kobs, when actin concentration is much larger than the myosin-S1 concentration ([A] ≫ [S1]). The absolute values of kobs calculated from the estimated KB and KT and obtained from fitting the observed stopped-flow transient at the range of Ca2+ concentrations were similar especially at lower Ca2+, whereas the calculated values overestimated the experimental fits at higher Ca2+ concentrations (Fig. 7A). The fractional 
kobsf obtained from the DLS, QN and SA methods showed a sigmoidal dependence on pCa (Fig. 7B) and differ little from the common sigmoidal shape even if KT fluctuate significantly for pCa< 6.0 (see Fig. 6B). Hill coefficients for kobs are also almost identical for the DSL, SA and QN methods, having values of 2.60, 2.61, and 2.71, respectively, and 50% of the change in kobs is achieved at virtually the same Ca2+ concentration (i.e. at pCa of 6.0). Interestingly, kobs is closely related to KB which is further exemplified by the mapping the spike in KB at pCa = 5.2 to the spike in kobs.</p><p>The kobs dependence of pCa derived from the best fits by any of the three-parameter estimation methods is similar to the kobs obtained from single exponential fits of experimental data (Boussouf et al., 2007a). Single exponential decay is observed at high Ca2+ concentrations and after a delay, tdelay, at lower Ca2+ concentrations. Fig. 8A shows an example of curve fitting of a single exponential at low Ca2+ concentration (pCa = 7.0). At low concentrations and after an initial delay, the single exponential (a straight line in semi-log plot), fits the data excellently until it reaches the noise floor at long times. These fits provide the single exponential rate constant kobs (s−1) that after appropriate normalization can be compared to the equivalent rate constant, kobs (or 
kobsf), calculated from estimates values of KB and KT (Eqs. (6) and (7)). Fitting the single exponential to the stopped-flow transient can be achieved by either fitting all available data or just a partial the range of data which excludes the initial delay before rapid decay and the data which reached the noise floor at long times. Fig. 8A demonstrates that more robust fits are obtained when the model is fitted to this partial range of data. Overall differences between the fits of all data vs. the partial range of data are small but significant (Fig. 8B), both fits are similar to data reported in (Boussouf et al., 2007a) and the delay time show the expected sigmoidal trend (Fig. 8C). However, the single exponential rate constant kobs underestimates the values of the equivalent rate constant at higher Ca2+ concentrations (Fig. 7A). This underestimation is likely due to the difficulty of separating the initial delay from the single exponential decay. This is also reflected in a lower Hill coefficient of 1.79 comparing to an average value of the coefficient of about 2.6 obtained from the estimated fractional equivalent rate constant 
kobsf (Fig. 7B).</p><p>We investigated the reason why the values of KB estimated by different methods are significantly different for pCa = 5.2 (Fig. 6A). We observed that the plot of the experimental data (Fig. 9) has a somewhat irregular shape during the first 60 ms that differs from the initial part of the sigmoidal curve observed at the other concentrations (inset in Fig. 9). For actin concentrations in excess of myosin-S1 ([A] ≫[S1]), this early-time irregularity is present at all concentrations but it is much smaller than at pCa = 5.2, only affecting the estimated KT values for pCa>6.0. It is interesting to note that despite quite different values of KB estimated by the DSL, QN and SA methods, the fitted curves are almost identical (Fig. 9), and the fit error is roughly the same (3.2603 × 10−5, 3.3183 × 10−5 and 3.3274 × 10−5, respectively). A possible reason for the differences observed in Kb is the relatively low value of the diagonal term of the resolution matrix associated to KB at pCa = 5.2 which is <0.2. Such low value indicates that the error of the fit is not sensitive to a relatively large variation of KB (Fig. 10A). Thus, the significantly different estimated values of KB are caused rather by different convergence behaviors of the three-parameter estimation methods than by the particular accuracy of the fit.</p><p>Overall, SA provided the most reliable estimates of KB and KT because of its ability to search for global minima in a broad parameter domain so that its estimates are virtually independent of the initial guess. This method provided stable convergence in all examined cases and the smallest fitting error. Interestingly, the DLS estimates agree with the estimates from SA in most of the cases considered in this study, which suggests that DLS can be used equally well as SA to fit data from stopped-flow fast transient experiments. In a few occasions DLS achieved even slightly better fits of the data but the differences in the error were negligible. In addition to an excellent fit, the DLS method provides additional information about the quality of the fit via the resolution matrix (RM). Fig. 10 shows the diagonal elements of the RM for KB and KT. Consistent with Table 2 and Fig. 5, these data indicate that KB is well resolved, especially for low calcium concentrations. The sharp decrease in the quality of the fit that is observed at high Ca2+ concentrations (pCa >6.0, RM of KB< 0.5) can be linked to the kinetics of the system (Fig. 10A). In fact, most of TmTn units are in either closed or open state at high Ca2+ concentrations (i.e. low pCa), thus further increase of KB causes a very small increase in the effective binding rate of S1 to regulated actin. In binding kinetic terms this means that at larger values of KB (>3) the overall rate binding of S1 to regulated actin is affected very little by opening of a few extra units (Fig. 10A). Consequently, large changes in KB values cause small shifts of the predicted transients which in extreme cases may be within the experimental uncertainty. At lower Ca2+ concentrations, however, KB is better resolved suggesting that the primary regulation mechanism is the transition between blocked and closed sates, i.e.by Kb.</p><p>In contrast to KB, KT is better resolved at higher Ca2+ concentrations (pCa>6.0, RM of KT>0.5) than at lower Ca2+ concentrations (Fig. 10B). This result suggests that the primary regulation mechanism at high calcium concentrations is the transition between closed and open states. This is because at higher Ca2+ concentrations the closed and open states are reasonably well occupied. However, the values of both KT and RM vary significantly from one Ca2+ concentration to the next, suggesting that these variations are more likely to be caused by imperfections in experimental methods and data acquisition than by inadequate parameter estimation. This argument is supported by Fig. 6B which displays the same pattern of seesaw variation of KT and essentially the same estimated Kt values with Ca2+ concentration regardless of the parameter estimation methods used. Consistent with this idea, the time courses of fractional fluorescence predicted by three-parameter estimation methods for the same Ca2+ concentration are very similar, although they significantly deviate from the observations (Fig. 9). Thus, the observed seesaw variation of Kt with Ca2+ concentration is not an artefact of the parameter estimation methods but rather a consequence of imperfections in experimentation to which Kt is very sensitive.</p><!><p>There are three desirable features of a good parameter estimation method: (1) ability to escape from local error minima; (2) uniqueness of the estimated parameters; and (3) fast convergence. The first of these features can be achieved with global search methods such as SA. However, as we showed above, some gradient methods, such as DLS, may fit the experimental data equally well. In the circumstances when gradient methods are equally effective in finding an absolute minimum of the estimation error, the gradient methods should be preferred because they converge faster, while at the same time providing information about the uniqueness of the solution via the resolution matrix.</p><p>Fig. 11 shows typical iteration histories of three MG model parameters (Kb, Kt and K2) and of the fit error at pCa = 4.6. This figure shows the estimation of the three parameters of the MG model because the estimation of all four free MG model parameters can only resolve well K1 having a value of about 0.25 μM−1, while other parameters are poorly resolved (see RM in Table 2 and Fig. 5). Thus, we fixed K1 to 0.25 μM−1 and analyzed the convergence of KB, Kt and K2 which can now be better resolved. Fig. 11 shows the convergence of KB, KT and K2, estimated by fitting the fluorescence transient for pCa = 4.6 and ε = 0.3. It is important to note that K2, KB and KT showed damped oscillations around their final convergence values. The interdependence between the estimated parameters is reflected in the RM (see Table 2, pCa = 4.6, three parameters and Fig. 5B), whose diagonal elements are appreciably smaller than 1 (KB ≈ 0.47, KT ≈ 0.58 and K2 ≈ 0.37), and whose off-diagonal elements corresponding to KT–K2 and in lesser degree Kb–Kt are relatively large. Note that off-diagonal of the RM elements indicate the degree of correlation between the parameters. For example, at pCa 4.6, KB can be well resolved if the number of parameters is reduced from three to two by fixing K2 to 170. In this case diagonal elements of RM have values >0.75, whereas offdiagonal elements are much smaller (~0.1, see Fig. 5C). The overall convergence can be improved by a better choice of an initial guess of the estimated parameters and by using a smaller value of the damping coefficient ε.</p><p>The effect of the chosen value of s on KB and the convergence and error of the three-parameter estimation (of KB, KT and K2) is shown in Fig. 12. As expected, the convergence is faster and the RM improves for lower values of ε. However, decreasing ε below 0.1 causes numerical instabilities that appear as ripples in the error function, and which could lead to divergence of the solutions. The same patterns are observed in the speed of convergence which rapidly increases for smaller values of ε. For example, the mean value of 
K¯B(r) converges very quickly for smaller values of ε (Fig. 13A) and the variations between neighboring iterations, 
ΔKB/K¯B(r) rapidly decrease (Fig. 13B). However, for the very small value of ε = 0.05, the approximate value of the converged 
K¯B(r) is achieved in just one or two iterations, but this value is not stable and the estimates kept oscillating around it without converging. In fact in most cases, when ε ≤ 0.1 smooth convergence cannot be achieved at all – unlike the case shown in Fig. 13. Fig. 14 shows that in this particular example the resolution matrix does not improve substantially when the damping parameter is decreased below 0.15. Note that the DLS method formally approaches the QN algorithm for very small values of ε, which in turn explains why the QN method is only effective when the initial guess is close to values of fully converged parameters. Thus, the use of the DLS method requires a systematic way of choosing initial guesses for the parameters to be estimated and a careful choice of ε. In our experience, the strategy that proved most successful in fitting the kinetic data is to start with a relatively large damping coefficient, say ε = 0.5, and gradually decrease ε while monitoring the convergence and the stability of the solution. The last estimated parameter values are the best initial guesses for the next parameter estimation with smaller ε values.</p><!><p>Estimating chemical kinetic rate constants from experimental data is complicated due to the nonlinear behavior of protein binding and because the data are prone to procedural imperfections and fluctuations. We developed algorithms to accomplish this goal with the McKillop–Geeves model of thin filament regulation in solution (McKillop and Geeves, 1993), and we suggested strategies of how to successfully extract the best set of values that enable to establish the Ca2+ dependence of the model parameters. Establishing relationships between the reaction constants (KB and KT) and the calcium concentration (pCa) is essential for a better understanding and further the development of quantitative models of thin filament regulation both, in solution and in living muscle fibers.</p><p>Three model parameter estimation methods, (i) Damped Least Squares (DLS), (ii) quasi-Newton (QN) which is related to DLS but simpler, and (iii) simulated annealing (SA), were compared by estimating the parameters of the McKillop-Geeves three-state rigid chain model from a set of stopped-flow data collected at different Ca2+ concentrations. All three models used different minimization methods to find best fits to a family of stopped-flow transients of S1 binding to regulated actin: (1) the DLS method is based on a gradient method for minimization of the error, (2) the QN method is a powerful optimization method close to DLS; and (3) SA is a global optimization method with an uphill and downhill search.</p><p>The MG model has four free parameters, KB, KT, K1, and K2, but the estimated parameters resulting from these methods have not provided any reasonable dependence of the model parameters on the Ca2+ concentration because the model parameters are interdependent and cannot be uniquely resolved (Table 2, and Figs. 5 and 10). Using evidence from the literature and from our own investigations we found that K1, and K2 only mildly depend on Ca2+ concentration and the trend is quite inconsistent. By setting reasonable values for K1 and K2, we found sigmoidal dependence of Kb on Ca2+ concentration (Hills coefficient of ~2.25) and a somewhat less prominent and scattered dependence of Kt on Ca2+. The parameters estimated by DLS, SA and QN methods have similar values for all of the studied Ca2+ concentrations (Fig. 6). The fractional 
kobsf which is the effective kinetic binding parameter, derived from model parameter KB and KT (Eq. (6)), also displays a sigmoidal dependence on Ca2+ concentration, and the difference between the parameters estimated by the different methods was minor (Fig. 6). Furthermore, the Hill coefficients of 
kobsf were virtually the same (2.6), and 
kobsf reaches 50% at pCa = 6.0.</p><p>Overall, SA appears to be the most reliable of the three estimation methods because it always converges and it is based on a global searching strategy. However, this method converges very slowly and requires a large number of model simulations. Thus, SA is limited for parameter estimation only if the model predictions can be calculated quickly so the converged solution can be achieved over reasonable computational time. In contrast, those models which explicitly include cooperativity between nearest neighbor TnTm units such as the Hill two-state model (Chen et al., 2001; Hill et al., 1980) or the MG cooperative model which can be only be solved by stochastic (i.e. Monte Carlo) simulations will require an over-extensive amount of computational time. Because these simulations are computationally intensive, only fast converging parameter estimation methods can be effectively used, i.e. the methods which converge in up to 25 iterations. One possibility is the use of the quasi-Newton method which converges rapidly, but its convergence is influenced frequently by the initial guess values and overall it requires a sizeable number of functional evaluations too. This method also imposes a large degree of uncertainty because choosing different initial guess values leads to different estimated model parameter values upon the incomplete convergence. The good fits shown in Figs. 3, 4, 6, 7 and 10, as well the parameter values shown in Table 1 are achieved by guessing initial values close to the values estimated by SA or DSL method; otherwise, the convergence of QN was compromised. These properties limit the application of the QN method to models based on stochastic simulations as well.</p><p>The parameter estimation via the sensitivity matrix and the Damped Least Squares (i.e. DLS) method has a number of advantages over the other two methods. Figs. 12 and 13 demonstrate that DLS can converge in less than 20 iterations, when an appropriate value of the damping parameter is used (i.e. ε <0.2), and reasonably well-chosen initial guesses are used. Similar to other iterative least squares methods, DLS can only find solutions that are linearly close to the initial guess and sometimes converge to a local minimum. These problems can be overcome by initially setting ε equal to a relatively large value (e.g. ε = 0.5) for a few iterations (say up to 20) and then progressively decreasing ε until convergence is reached. The initial large value of ε will allow a wider search within the parameter space but it may not converge. Typically, the values of each of the estimated parameters as well as the error fluctuate between some larger and smaller value between two consecutive iterations (Fig. 11). In most cases the estimated parameters will not converge in 20 iterations, but the last few iterations provide sufficient information to choose a new guess which is much closer to the converged solution. The criterion to select new initial parameters should include values with the lowest error. The subsequent estimation with the new guess and ε < 0.2, yields rapid convergence. This methodology was successfully applied to the estimation of the MG model parameters. It provided almost the same parameter estimation values as the global search method SA over all Ca2+ concentrations, but in fewer iterations (i.e. <20). Note that the strategy presented above (starting with lager ε and reducing it as approaching a minimum) is a standard optimization procedure for a Levenberg–Marquardt method (Kecman, 2001; Nocedal and Wright, 2006).</p><p>The uniqueness of the estimated parameters is a very important point to consider when dealing with parameter estimation methods. It is well-known that multiple solutions may exist with equal error margins (Tang et al., 2005; Vajda and Rabitz, 1988; Zsely et al., 2004). Various combinations of parameters can fit the experimental data equally well within the experimental uncertainty, or some of the parameters may be interdependent. Among those multiple solutions, the "global" solution might be attractive by offering the smallest residual error. Because we are dealing here with a nonlinear inversion problem, it is desirable to examine the global properties of the fit error in order to be confident that the most accurate set of parameters is estimated. In our case, SA is a convenient method for exploration of the global parameter space searching for the absolute minimum error. However, this global residual error minimum is not guaranteed to be better than other local minima estimated by other methods such as DLS (Tang et al., 2005). This observation is indeed true in the current instance (see Table 1 and Fig. 5), where KB–pCa relationships obtained by SA and DLS differ minimally. They only differ appreciably at high Ca2+ concentrations where the error landscape is relatively flat and a wide range of estimated KB values produces model predictions with virtually no change in error. Thus, in a noisy and relatively flat error landscape the meaning of the absolute minimum is somewhat diluted and in this situations there is no a significant benefit in using global search methods, such as SA. In this case, the RM obtained from DLS can provide additional information on how to obtain the best solution with a minimum number of iterations by adjusting the appropriate value of ε and monitoring the variation of estimated parameters and the RM (Table 3, Figs. 12 and 14). In general the smaller the damping parameter ε, the faster is the convergence and the model parameters are better resolved. However, too small ε may lead to numerical instabilities and divergence of the parameter estimation (Figs. 12 and 13). As pointed out by Menke (1989) choosing an appropriate value for ε is not a straight forward process and typically ε is set by trial and error in usually decreasing sequences of values.</p><p>In conclusion, if the complex error landscape has multiple minima, then the solution can be non-unique and a priori information must be added to resolve the indeterminacy. Therefore, there is no guarantee that the DLS or any other iterative technique will converge to the solution, even if it is known that the unique solution exists, because it is necessary to explore the global error surface to ensure that a global minimum is reached. Overall the SA method is the most reliable method for the estimation of parameters in complex chemical systems. However, this method is in some cases impractical because it requires an extremely large number of iterations to find absolute minimum of the error. Thus, other methods, such as DLS, or sophisticated combinations of methods should be used. Nevertheless, DLS parameter estimation of the MG model parameters provided almost identical parameter estimates as the SA method and supplies information about the interdependence of the estimated parameters. Furthermore, the main advantage of the DLS vs. other examined methods is its fast convergence which is essential for the solvability of parameters in complex reaction systems which can only by simulated by complex computational algorithms. It should be noted that these conclusions are specific to the single model of myosin binding kinetics to regulated actin. The relative advantage of DSL observed in this study should be further evaluated for other models of mechanochemical systems. It would be especially important to establish the effectiveness of parameter estimation methods when applied to Monte Carlo simulations such as Hill's model of thin filament regulation (Chen et al., 2001; Hill et al., 1980) or Continuous Flexible TnTm Chain model (Smith and Geeves, 2003; Smith et al., 2003). The ultimate parameter estimation method that employs these stochastic simulations would require: (1) fast convergence because the Monte Carlo simulations require extensive computational time for each simulation and (2) sufficiently robust algorithm which can uniquely resolve complex issues related to interference between noise in experimental data and inherent stochastic noise in simulation data.</p>