A novel kinetic model to demonstrate the independent effects of ATP and ADP/Pi concentrations on sarcomere function

Understanding muscle contraction mechanisms is a standing challenge, and one of the approaches has been to create models of the sarcomere–the basic contractile unit of striated muscle. While these models have been successful in elucidating many aspects of muscle contraction, they fall short in explaining the energetics of functional phenomena, such as rigor, and in particular, their dependence on the concentrations of the biomolecules involved in the cross-bridge cycle. Our hypothesis posits that the stochastic time delay between ATP adsorption and ADP/Pi release in the cross-bridge cycle necessitates a modeling approach where the rates of these two reaction steps are controlled by two independent parts of the total free energy change of the hydrolysis reaction. To test this hypothesis, we built a two-filament, stochastic-mechanical half-sarcomere model that separates the energetic roles of ATP and ADP/Pi in the cross-bridge cycle’s free energy landscape. Our results clearly demonstrate that there is a nontrivial dependence of the cross-bridge cycle’s kinetics on the independent concentrations of ATP, ADP, and Pi. The simplicity of the proposed model allows for analytical solutions of the more basic systems, which provide novel insight into the dominant mechanisms driving some of the experimentally observed contractile phenomena.

bound to a 1 ) and myosin node #2 is bound to actin node #2 (m 2 bound to a 2 ).

Z-line node:
k The following example of vectors P , vector V , and matrix K represent the linear equations for a 4-node system where there are no cross-bridges binding the actin and myosin filaments.Force output is determined by the external spring element bound to the M-line of the half-sarcomere.The titin element binds the Z-line to the first myosin node, and its rest length is denoted by γ 0 .For simplicity, we assume titin exerts a passive force when the half-sarcomere is stretched beyond its rest length, but does not exert a compressive force when the half-sarcomere is shortened.A spring constant of 10 pN/nm was chosen for titin because only very small stretching deformations are expected in this model [6].If the model is to be utilized with large stretches of the sarcomere from its rest length, then a nonlinear spring constant for k T should be used instead [6].
July 18, 2024 Supporting Information S2/S15 Equation S7 can be re-written for the situation where the cross-bridges are engaged between nodes a 1 and m 1 and between nodes a 2 and m 2 .Both cross-bridges are in their bound, pre-power stroke conformation.Force output is determined by the external spring element bound to the M-line of the half-sarcomere: The following equations describe a binding scheme identical to the previous set of equations.However, node m 1 is now in state 3, its post-power-stroke, high force bearing July 18, 2024 Supporting Information S3/S15 state.This state is characterized by a conformational change, which changes the rest length of the m 1 myosin head.This change manifests as an adjustment in the b 0 terms in the V matrix for the nodes involved in the state 3 cross-bridge i.e. a 1 and m 1 , changing these terms to include (b 0 − d ps ) instead, where d ps is the length of the power stroke.The power stroke in this model generates force in the positive x-direction:

C. Implementation of Kinetics
The stochastic nature of the three-state cross-bridge cycle used in this study was modeled using a Monte Carlo algorithm, more specifically, we implemented an algorithm similar to Gillespie [7].To describe our algorithm, we first have to define all rate constants, and for that we have to define all relevant free energies (see also Remaining variables are defined in Table 1.
Rate constants of transitions are detailed in Equations 6-13.Consistent with basic principles, the overall reaction quotient, depends only on the combination of concentrations [ATP]/[ADP][Pi] and hydrolysis free energy ∆G hyd , while individual rate constants depend on concentrations beyond this usual "mass action" combination.
Having defined rate constants, we advance time in (uneven!) steps δt, and the probability of the i → j transition during one such time step should be P ij = k ij δt.We have to realize that any one such transition in the many-myosin system affects rates of all possible subsequent transitions.This is because of the elastic coupling between myosin heads, expressed in the transition rates' dependencies on the potential energies and coordinates, x.Therefore, to faithfully represent the correct Poisson statistics of the molecular transitions, it is imperative to choose time step δt so small that only one transition can occur during that time step.To achieve this while avoiding the computationally prohibitive slowness of the algorithm, we vary time steps as follows (see The main idea is to choose the time step such as to limit the chance of an "uneventful" time step, when no transition occurs.We, therefore, choose (arbitrarily) the no-transition probability as a parameter, P no transition ; we chose P no transition = 0.1.In other words, some transition should occur with probability 1 − P no transition = 0.9.With this in mind, we sum up all the rates of all possible transitions from a current state and Simulation-representative probability bar of the system above.The size of the probability bins are proportional to their true values as seen in a system where cross-bridges are perfectly aligned with an actin node i.e. there is no deformation of the cross-bridge springs.Note that the probability of no transition occurring is 10% for every time step, regardless of the rate constants' values.
find the next time step as follows: Here m total is the number of myosin heads in the system; k n ij and k n ik are, respectively, forward and backward rates of myosin number n, which presently is in the state i.
July 18, 2024 Supporting Information S6/S15 Then, the probability of i → j transition for myosin n is expressed as Knowing this, we separate the interval (0, 1) into 2 * m total + 1: two intervals for each myosin with lengths P n ij and P n ik , plus one interval of length P no transition .We then generate a random number p * from a uniform distribution between 0 and 1 using MATLAB's internal rand() function (reshuffling random number generator at the start of each simulation using the rng('shuffle') function).Depending on the value of this random number, we choose the transition to perform.
Note that although p * is uniformly distributed between 0 and 1, it does not represent the probability distribution of state transitions.The latter is dictated by the rate constants and, therefore, changes with internal deformations of the sarcomere model.
Probabilities of state transitions occurring are calculated as the product of the rate constant and time step (Equation S17).After state transitions are defined by the Monte Carlo algorithm, we assumed mechanical equilibrium is achieved rapidly relative to the time taken to transition between states.The position of each node within the lattice was calculated by solving a detailed force balance at each node (Equation S10).The formulation of the model is flexible enough to allow for the addition of more actin binding sites (actin nodes) and cross-bridges (myosin nodes).

D. Parameter Exploration
Parameter exploration was performed on this sarcomere model to determine an appropriate set of values for the rate constants.The output of the model to which the parameters were fit was the myosin duty ratio, which in striated muscle is on the order of 0.05 [12][13][14][15].The myosin duty ratio is defined as the proportion of the cross-bridge cycle that the motor domain of the myosin is strongly bound to the actin filament [12][13][14].For this parameter exploration, the duty ratio was calculated can be found in Table 1.These potential parameter values were varied in a systematic manner to explore the full range of possibilities for each parameter.To determine the optimal parameter values, the model was run for each set of three parameters and the resulting myosin duty ratio was calculated (Fig B -D).The set of parameter values that resulted in the best fit to the duty ratio was selected as the optimal set.The selected parameters for the simulations in this study were well within the blue region of the plots in Fig B -D, meaning that the results of the model are not sensitive to slight parameter perturbations.
The task of this study was not to find optimized sets of model parameters that will apply to all scenarios, but rather to demonstrate how outputs of the model change in response to different geometries and metabolite conditions under a single set of parameters.Thus, parameter values that were chosen for the model were on the order or within physiological ranges with reasonable accuracy.A sensitivity analysis was performed where small perturbations in all of the parameter values were analyzed.The model was found to be robust to such perturbations for the outputs of duty ratio, ATP consumption rate, and average force output.Furthermore, because literature reported values of duty ratio exist within a range (3-5%), the model's parameters can be adjusted in ways to achieve different duty ratios or geometries associated with different muscle types (e.g.fast, slow, cardiac), myosin classes, and species [3,13,14].Similarly, mechanical parameters can also be adjusted in order to model sarcomeres from different physiology.For example, thick filament stiffness varies between different species, and this parameter can be adjusted to match these differences between species [16].In another example, equilibrium constants between cross-bridge states, and thus free energy changes between states, also vary between myosin isoforms [14,[17][18][19][20].  [3,21,22,[25][26][27][28][29].Additionally, literature reports of how many myosin actually participate in contraction (50-60% of myosin are in the SRX state and do not participate in contractility) was also taken into consideration [30].Lastly, some of the experiments from which estimates of ATP consumption were taken were conducted at lower than physiological temperatures and non-maximal calcium activation.Increased temperature and activation of thin filament binding sites would further increase ATPase rates in these systems.The goal of this study concerning ATP consumption was not to optimize it for all situations, but to observe how this activity varies under different conditions.

Fig 3 and
Equations 2-5 in the main text).State 1 was assigned a reference energy of 0 (Equation 2).The remaining free energy July 18, 2024 Supporting Information S4/S15 profiles of state 2 and state 3 are parabolic in relation to the distortion (x m − x a − b 0 ) of the cross-bridge in those states, where x m represents the x-coordinate of the myosin node, x a represents the x-coordinate of the associated actin node (nearest, unoccupied node), and b 0 represents the rest length of a cross-bridge.∆G * hyd is the standard free energy of ATP hydrolysis.The power stroke distance is represented by d ps , k B is the Boltzmann constant, T is absolute temperature, assumed to be room temperature at 300K.Under standard conditions, [ATP] * = [ADP] * = [Pi] * = 1 M.The normal concentrations for each of these molecules used in this study are: [ATP] = 5 mM, [ADP] = 0.03 mM, and [Pi] = 3 mM [8-11].
analytically from the equilibrium constants of each individual cross-bridge transition step (Equations 16, 19-21).These analytical values were found to match those achieved via simulation (Fig 5).Value ranges for the rate constants were taken from previous models and physiological values from literature.Literature values for these parameters July 18 Fig B. (Part1) Parameter exploration for constants in rate equations.Each constant was evaluated within its relevant range against two other constants.The model output used to fit the parameters was the myosin duty ratio, which in skeletal muscle ranges from 3-5% based on previous studies.

1 )
Fig C. (Part 2) Parameter exploration for constants in rate equations.Each constant was evaluated within its relevant range against two other constants.The model output used to fit the parameters was the myosin duty ratio, which in skeletal muscle ranges from 3-5% based on previous studies.

m 1 st. 1 to st. 2 m 1 st. 1 to st. 3 no trans. m 2 st. 1 to st. 3 m 2 st. 1 to st. 2 0 1 p* m 1 st. 1 m 2 st. 1 (C) Schematic of the Probability Bar for Myosin 1 in state 2 + Myosin 2 in state 1 m 1 st. 2 to st. 3 m 1 st. 2 to st. 1 no trans. m 2 st. 1 to st. 3 m 2 st. 1 to st. 2 0 1
Fig A. Probability bar depicting implementation kinetics using Monte Carlo algorithm.(A)Schematic of how the probability bar may appear for a 2-myosin system where both myosin are in state 1.Half-sarcomere schematic depicts the physical representation of this system.(B) A random number p * is pulled from a uniform distribution between (0, 1).In this example, p * falls into the first bin of the probability bar, indicating a forward transition of m 1 from state 1 to state 2. (C) Schematic of the probability bar following this transition, where m 1 is now in state 2 and m 2 is still in state 1. (D)