3  Nanoscale: modeling molecular motors

In Chapter 2, we presented a modeling framework of the cardiac tissue relying on the additive decomposition of the total stress into a passive and an active contribution. The active stress in the direction of the fiber was written \[ \underline{\underline{\Sigma}}_{\raise 0.3em {\,a}}= \frac{T_{\text{fib}}}{(1+2\underline{\tau}\cdot\underline{\underline{e}}\cdot\underline{\tau})^{\frac{1}{2}}}\underline{\tau}\otimes\underline{\tau},\quad\text{with } T_{\text{fib}}= \rho_{\text{m}}\tau_{c}, \] where \(\rho_{\text{m}}\) is the number of myosin molecular motors in a layer of thickness \(\ell_{\text{fib}}\) per unit cross-sectional area, and \(\tau_{c}\) is the force produced by a single motor. This chapter summarizes our contribution to the modeling of this force.

Associated references

The results summarized here are from the work of Caruel, Moireau, and Chapelle (2019), Kimmig and Caruel (2020), Chaintron, Caruel, and Kimmig (2023), and Chaintron et al. (2023).

We recall that the large majority of existing heart models are based on a direct rescaling of the average force generated by an infinitely large population of independent molecular motors, which is equivalent to rescaling a single representative motor.

3.1 Background: the myosin molecular motor

3.1.1 Myosin structure

Figure 3.1: Schematic representation of a myosin II heavy chain. A myosin II protein consists of two heavy chains coiled together by their tails. Each chain contains a motor domain that interacts with actin and that contains the active site of ATP hydrolysis. The force is genarated upon conformational changes of the motor domain that gets amplified by the lever arm structure. Figure taken from (Robert-Paganin et al. 2020)

The myosin II protein responsible for muscle contraction consists of two heavy chains containing about 2000 amino-acids each. The two chains are coiled together by their tail regions while their head regions remain independent. A schematic representation one of the chains is shown in figure Figure 3.1.

Figure 3.2: Ribon representation of a molecular structure with a myosin VI molecule bound to actin ontained using Cryogenic Electron Microscopy (Cryo-EM). Figure produced by Manevy (2023) based on the structure 6BNQ of the protein databank.

The motor domain is the part of the protein that interacts directly with the actin filament during the contraction. It is connected to a coiled-coiled part of the protein called the lever arm. Its christallographic structure has been resolved in the 90s (Rayment et al. 1993; Dominguez et al. 1998), see also the historical review presented in Robert-Paganin et al. (2020) for other relevant references. A structure of the motor domain of a myosin head (without the lever arm) bound to actin shown in Figure 3.2. The motor domain itsel is decomposed into several subdomains that surround the active site, where the ATP hydrolysis occurs, see Figure 3.3.

Figure 3.3: Repsentation of subdomains of myosin motor domain and lever arm surrounding the active site where the ATP hydrolysis products, ADP and Pi are located. Figure taken from (Robert-Paganin et al. 2020)
Figure 3.4: Structure of the pre- and post conformation of myosin II. The superimposed motor domains are highlighed in yellow. The pre-power stroke conformation has the upward lever-arm (pink, blue and cyan). The total displacement between the ends of the two lever arm configuration is about 12 nm. Figure taken from (Dominguez et al. 1998)

An important aspect of the force generation mechanism is the ability of the myosin motor domain to undergo a conformational change that is amplifies by the lever arm in order to induce the relative displacement of the myosin and actin filaments. This conformational change is called the power stroke (or working stroke). The molecular structures characterizing the power stroke has been resolved by Dominguez et al. (1998). A representation of the pre- and post-power stroke structures is shown in Figure 3.4.

3.1.2 Structural aspects of the myosin-actin interaction cycle

In our Introduction, we presented a simplified version of the Lymn and Taylor (1971) cycle, see Figure 1.3. This cycle has been refined since then based on numerous structural studies that revealed the structures corresponding to various intermediate steps. As an example, the cycle proposed by Houdusse and Sweeney (2016) is shown in Figure 3.5.

Figure 3.5: Representation of the actin-myosin interaction cycle emphasizing the state of the motor domain and the position of the lever arm. Figure taken from (Houdusse and Sweeney 2016)

While structural studies undeniably offer valuable insights into the molecular processe of force generation, they do not capture the dynamics of the transitions between the states. However, these structures can serve as initial confitions of Molecular Dynamics (MD) Simulations that can be performed to compute the associated free energy landscapes and short timescale dynamics (Fischer et al. 2005; Baker and Voth 2013). Different levels of description are available: all atoms or with different degrees of coarse-graining (Tozzini 2005). These costly simulations, even the most coarse-grained ones, can cover timescales up to ~1 μs , which is ineffective to simulate the conformational changes involved in force production which last more than 1 ms (Huxley and Simmons 1971). Moreover, current MD studies focus on single transitions of the actin-myosin interaction cycle (see Figure 3.5). Simulating the entire ATPase cycle for even a single motor therefore remains a significant challenge.

These simulations however provide valuable information on the effective mechanical properties of the proteins and help to determine the collective variables that can be used to caracterize the different transitions within the cycle. Starting from an all atom dynamics in a space of very high dimension, the strategy is to project this dynamics on lower-dimension manifold parametrized by a set of collective variables to obtain a reduced model.

The highest degree of simplification is when the reduced dynamics becomes a jump process between discrete states, as depicted in Figure 3.5. This description thus views the actin-myosin interaction cycle as a sequence of chemical reactions whose rates determine the dynamics of the system. Most of the currently used molecular motors models are based on this idea, see Section 3.2 for examples.

3.1.3 Single molecule force spectroscopy experiments

Nanomanipulation techniques such as optical tweezers, Atomic Force Microscopy or Scanning probe microscopy, have been developed since the 1990s to study mechanostransduction phenomena of the living cells and in particular force generation by molecular motors (Spudich et al. 2011). They have been particularly successful at determining the functioning of processive molecular motors, i.e. motors that operate individually (Kinesin, Myosin VI etc., see (Block, Goldstein, and Schnapp 1990; Svoboda et al. 1993; Rock et al. 2001; Ménétrey et al. 2005)) or slow muscles myosins (Veigel et al. 2003). Unlike muscle Myosin II, processive molecular motors stay permanently in contact with their track thanks to their pair of heads actin like coordinated legs. Non processive motors like muscle Myosin II spend only a small fraction of their cycle attached to their track. Studying these short-lived attached configuration requires sub-millisecond resolution of the loading device, which is a signigicant technical limitation (Finer, Simmons, and Spudich 1994; Molloy et al. 1995; Veigel et al. 2003).

The development ultrafast force clamp single molecule spectroscopy developed over the past decade, (see in particular Capitanio et al. 2012), alows now to detect events occuring at timescales comparable to that of the power stroke (Woody et al. 2019), which partially remove the previous limitations. Before that however, we can mention the work of Veigel et al. (2003) on the slower smooth muscle myosin who showed that the myosin detachment rate is load-dependent: it is larger when the molecule is in compression than when in tension.

A contemporary promissing application of single-molecule biophysical experiments is the in-vitro testing of mutated proteins and specific modulators used as treatments of muscular diseases. A recent example can be found in (Woody et al. 2018) where a dual-laser optical trapping setup was used to study the effect of Omecamtiv Mecarbil, a small positive cardiac inotropic agent that increases cardiac performance in acute heart failure situations (Malik et al. 2011).

Since our work is focused on the modeling of non-processive motors, the systematic comparison of the model predictions with single modelule experimental results was not pursued, though it is not excluded for future model calibration procedures.

3.2 Existing models

3.2.1 Chemical-mechanical models

Associated references

The following is a presentation of the most widely used class of actin-myosin interaction models. While the orginal formulation dates back to the 1950s and 1970s, more recent presentations with additional mathematical and numerical developments can be found in (Kimmig, Chapelle, and Moireau 2019; Chaintron, Caruel, and Kimmig 2023). The following paragraphs are adapted from these references.

Figure 3.6: The original two-state Huxley–Hill model (Huxley 1957). The model is parametrized by the distance \(s(t)\) separating the myosin anchor from its nearest binding site. The distance between two consecutive sites is denoted by \(d\). After attachment, the formed cross-bridge is represented as a spring of length \(s(t)\). Adapted from (Chaintron, Caruel, and Kimmig 2023)

The typical reprentation of the actin-myosin dynamics is a discrete jump process between states endowed with mechanical degrees of freedom. The first model of this kind was proposed by A.F. Huxley in 1957 (Huxley 1957). It has then been formalized and generalized by T.L. Hill and co-workers in 1974 (Hill 1974). Its formulation is illustrated in Figure 3.6 1. The model is parametrized by the distance \(s(t)\) separating the myosin anchor from its nearest binding site. The distance between two consecutive sites is denoted by \(d\). A fundamental hypothesis is that a given head can interact only with the nearest binding site. A myosin head exists in two states: detached or attached to actin, whith attachment and detachment events being modeled as a “chemical” jump process with rates that depend on the distance to the nearest binding site \(s\). The Huxley’ 57 model sets the foundation for all subsequent models that use this so-called chemical-mechanical approach.

1 In (Chaintron, Caruel, and Kimmig 2023), we discuss the well-posedness conditions of the Huxley-Hill model and propose an alternative formulation based on the use of Poisson random measures. We determine conditions that ensure the model’s compatibility with the thermodynamic principles.

The collective variable that lumps the high dimensional space of protein configurations is here the discrete variable \(\alpha_{t}\in\{0,1\}\). This variable characterizes the attachment state of the myosin head, taking the value \(\alpha_{t}=1\) if the head is attached, and the value \(\alpha_{t}=0\) if the head is detached. After attachment, the formed cross-bridge is represented as a spring of length \(s(t)\), see Figure 3.6.

The model assumes a large population of independent myosin proteins having their closest attachment site located at \(s \in (-d/2,+d/2)\). The function \((t,s) \mapsto P_1(t,s)\), represents the fraction of this population that is attached at time \(t\). Similarly, the function \((t,s) \mapsto P_0(t,s) := 1-P_{1}(s,t)\) represents the fraction of this population that is detached at time \(t\).

Figure 3.7: Stochastic jump process dynamics and energetics. (a) illustration of the transition rates between the attached and the detached states. (b) illustration of the first principle. Adapted from (Chaintron, Caruel, and Kimmig 2023).

The stochastic jumps dynamcics associated with \(\alpha_{t}\) involves four transition rates (see Figure 3.7 (a)):

  • \(k_{0\rightarrow 1}(s)\) associated with the jump \(0 \rightarrow 1\) (direct attachement);
  • \(k^{\mathrm{rev}}_{0\rightarrow 1}(s)\) asociated witht the jump \(1 \rightarrow 0\) (reverse attachement);
  • \(k_{1\rightarrow 0}(s)\) associated with the jump \(1\rightarrow 0\) (direct detachment);
  • \(k^{\mathrm{rev}}_{1\rightarrow 0}(s)\) associated with the jump \(0\rightarrow 1\) (reverse detachment).

These rates can be aggredated by defining \(f(s) = k_{0\rightarrow 1}(s)+k^{\mathrm{rev}}_{1\rightarrow 0}(s)\) and \(g(s) = k_{1\rightarrow 0}(s) + k^{\mathrm{rev}}_{0\rightarrow 1}(s)\). The fraction of attached myosin heads then verifies the following PDE, for \(t > 0, \, s \in \left(-\tfrac{d}{2},+\tfrac{d}{2}\right)\), \[ \left\{ \begin{aligned} & \partial_t P_{1}(t,s) + \dot{x}_{c}(t) \partial_s P_{1}(t,s) =-g(s) P_{1}(t,s) + f(s) \left[1 - P_{1}(t,s)\right], \\ & P_{1}(t,\pm d/2) = 0,\; t \geq 0,\\ &P_{1}(0, s) = P^\mathrm{ini}_{1}(s),\; s \in \left[-\tfrac{d}{2},+\tfrac{d}{2}\right], % & P_{0}(s,t) =1 - P_{1}(s,t)&. \end{aligned} \right. \tag{3.1}\] where \(\dot{x}_{c}\) is the relative sliding velocity between the two filaments, see Figure 3.6.

The tension generated by the attached motors at a gven sliding velocity \(\dot{x}_{c}\) is then defined as \[ \tau_c(t) := \frac{1}{d} \int_{-d/2}^{+d/2} \partial_s w_1(s) P_{1}(t,s) \mathrm{d}s, \] where \(w_{1}\) is the energy of the protein in the attached state. If the attached protein is viewed as an elastic spring, then \(w_{1}(s) = (\kappa/2)(s-s_{0})^{2}\), where \(\kappa\) is the stiffness and \(s_{0}\) the reference length, see Figure 3.7 (b). This tension is finally rescaled to obtain the macroscopic active stress \(T_{c}\) that appears the macroscopic beharior law, see Section 2.1.3 and Equation 2.2 in particular.

Overall, the thermodynamic properties of the models are defined by the energy landscapes \(w_{\alpha}\) and by the transition rates. There is no a priori guarantee that for an arbitrary choice of these paramters, the system would be compatible with the principle of thermodynamics These energetic aspects of the chemical mechanical models have been clarified in the wotk of Hill (1974) and Hill (1976). In particular, it was shown that the system satisfies the second principle of thermodynamics if for each transition between two states \(A\) and \(B\) of the cycle, the associated transition rates \(K_{A\to B}\) and \(K_{B\to A}\) verify the detailed balance condition \[ \frac{K_{A\to B}}{K_{B\to A}} = \exp\left[-\frac{w_{B}-w_{B}}{k_{\text{B}}T}\right], \] where \(w_{A}\) and \(w_{B}\) are the energies of states \(A\) and \(B\), respectively.

In the case of the two state model illustrated in Figure 3.7 (a), this property reads \[ \begin{aligned} \frac{K_{0\rightarrow 1}(s)}{K^{\mathrm{rev}}_{0\rightarrow 1}(s)} &= \exp\left[-\frac{w_{1}(s) - w_{1}(s)}{k_{\text{B}}T}\right],\\ \frac{K_{1\rightarrow 0}(s)}{K^{\mathrm{rev}}_{1\rightarrow 0}(s)} &= \exp\left[-\frac{(w_{0}(s)-\mu_{T})-w_{1}(s)}{k_{\text{B}}T}\right] \end{aligned} \tag{3.2}\]

Notice that, in the second equation, the energy of the detached state is \(w_{0}-\mu_{T}\) and not \(w_{0}\). The difference is the free energy \(\mu_{T}\) harvested from the hydrolysis of ATP. If \(\mu_{T}=0\), then it can be shown that the system reaches thermal equilibrium with \(\frac{P_{1}(s)}{1-P_{1}(s)}=\exp[-\frac{w_{1}(s)-w_{0}(s)}{k_{\text{B}}T}]\), and that there is no net flux in the cycle.

A typical cycle is shown by the thick line in Figure 3.7 (b), where we illustrate for each sequence of that cycle, the associated first principle quantity: work (\(\mathcal{W}\)), heat (\(\mathcal{Q}\)) and energy input (\(\dot{\mathcal{E}}\)). Having \(\mu_{t}>0\) ensures that, on average the overall free energy decays at each cycle, which maintain non-equilibrium. This point is discussed in a general context by Jülicher, Ajdari, and Prost (1997), see also (Jülicher and Prost 1995; Wang and Oster 2002).

The class of chemical-mechanical models was intensely studied in the late 90’s and up to the late 2010’s, in particular in relation with the development of biomimetic in vitro experiments (Holzbaur and Goldman 2010; Pollard 2018). These works shed light on a variety of dynamical behaviors akin to spontaneous oscillations, emerging particularly when a cluster of motors is coupled to a common elastic substrate. We refer to the work of Jülicher and Prost (1995), Jülicher, Ajdari, and Prost (1997), Duke (1999), Plaçais et al. (2009), Guérin et al. (2010) and Guérin, Prost, and Joanny (2011) for further details.

Today, the vast majority of existing models pertains to this category and differ from the original Huxley-Hill two state formulation essentially by the number of chemical states involved and by the nature of the spring (linear or nonlinear) (Huxley and Simmons 1971; Piazzesi and Lombardi 1995; Smith and Geeves 1995a, 1995b; Duke 1999; Smith et al. 2008; Smith and Mijailovich 2008; Månsson 2010, 2016, 2020; Caremani et al. 2015; Pertici et al. 2018).

In summary, the chemical-mechanical models projects the high dimentional molecular dynamics on a discrete space of states. The originality of the approach vis a vis a purely chemical system is the association of each state with a mechanical degree of freedom (here denoted as \(s(t)\)). The resuling energy landcape is thereby coupled with the relative position of the two filaments.

3.2.2 Stochastic ratchet models

An alternative to the Huxley-Hill-type models is to consider a less drastic dimensional reduction of the molecular motors’ molecular dynamics, by keeping a description of the “internal motion” of the proteins. Using such approach, the state of the motor is parametrized by a vector of continuous stochastic (collective) variables \(\boldsymbol{X}_{t}\) evolving in a low dimensional space. The dynamics of \(\boldsymbol{X}_{t}\) follows the gradient of an energy landscape wich represent the projection of the all-atom potential on the lower dimension space of collective variables.

A prototype model this kind was proposed by Magnasco (1993). The molecular motor is a particle \(X_{t}\in\mathbb{R}\) evolving on a unidimensional periodic energy landscape according to langevin stochastic differential equation

\[ \mathrm{d}X_{t} = -[\phi'(X_{t}) + f(t) + q]\mathrm{d}t + \sqrt{2D}\mathrm{d}B_{t}. \tag{3.3}\]

In Equation 3.3, \(\phi\) is a non-symetric periodic potential mimicking the interaction with the actin track, \(q\) is the external load, \(D\) is the diffusion coefficient and \(\mathrm{d}B_{t}\) is a brownian motion increment. The term \(f(t)\) represent a so-called colored component of the noise, as opposed to the pure white noise represented by the term \(\mathrm{d}B_{t}\). In the case of Magnasco (1993), it is represented by a deterministic piecewise constant force of the form \(f(t) = A(-1)^{n(t)}\) with \(n(t)=\lfloor 2t/\tau\rfloor\). The average of this force over a period \(\tau\) is equal to 0, implying that \(f\) has overall no net effect, but combined with the uncorrelated brownian motion and the non-symetric periodic landscape, it can effectively generate a directional motion in the opposite direction of the external force. This is a way to represent the chemical energy input from the ATP hydrolysis as a purely mechanical element.

3.3 Challenges

At the scale of the molecular motor several challenges can be stated.

  1. Not all the crystallographic structures associated with the different steps of the Lymn–Taylor cycle are resolved, and MD simulations only provides very short time dynamics, so one cannot use this technique to simulate the entire interaction cycle. The challenge here is to formulate new methodologies allowing to use the data generated by these simulations to formulate simplified models of the actin-myosin interaction, like for instance the chemical-mechanical models.

  2. A salient feature of the existing chemical-mechanical models is the use of functions to parametrize the transition rates between the states, see \(f\) and \(g\) in Equation 3.1. These rates are difficult to calibrate in an unequivocal way: they represent a virtually infinite set of parameters. Currently, the data available to constrain their calibration come mostly from experiments dealing with many motors in interaction, different animal species and heterogenous experimental conditions. Other data can be extracted from single molecule experiments. One could leverage the results of the molecular dynamics simulations on single proteins to complement experiments involving many motors.

  3. The calibration is greatly facilitated by the mean-field assumption mentioned above, since it implies a rather direct coupling between the experimental observables and the single motor characterisitcs. This fundamental assumption can be questionned.

3.4 Contributions

3.4.1 The Huxley-Hill formulation revisited

Associated reference

This section is an abstract of the model and results presented in (Chaintron, Caruel, and Kimmig 2023).

Figure 3.8: New formulation of the Huxley-Hill model proposed by Chaintron, Caruel, and Kimmig (2023). In this new formulation, called the h-model , the dynamics of the position of the detached head \(X_t\) relative to the same anchor point is described. The position of the nearest binding site, denoted by \(h_{t}\) is now defined with respect to this position \(X_{t}\) and not anymore with respect to the anchor of the myosin head in the thick filament, see Figure 3.6 for comparison.

Recently Chaintron, Caruel, and Kimmig (2023) a published a slightly modified formulation of the Huxley-Hill model, see Figure 3.8. In this new formulation, called the h-model, the state of the myosin motor is characterized by three stochastic variables:

  • \(X_{t}\), the position of the head with respect to the myosin enchor on the thick filament,
  • \(\alpha_{t}\), the attachment state (\(\alpha_{t}=0\) is the motor is detached and \(\alpha_{t}=1\) if it is attached),
  • \(h_{t}\), the position of the nearest actin site with respect to \(X_{t}\).

The main difference with the original formulation is that the position of the nearest biding \(h_{t}\) site is not anymore defined with respect to the anchor on the thick filament, see Figure 3.6. The two models are linked formally by the relation \(s_{t} = X_{t} + h_{t}\).2 Upon attachment, the variable \(h_{t}\) jumps to 0 so that \(X_{t}=s_{t}\) after attachment.

2 Note that here \(s(t)\) is the distance between the anchor on the myosin filament and the site closest to \(X_t\), which is different from the Huxley-Hill model.

The consequence of this seemingly innocent change is that the myosin head can now attach to any site on the real axis without having to resolve to an explicitly multi-site formulation like for instance in (Kimmig, Chapelle, and Moireau 2019; Månsson 2019). Using a multi-site model comes at the price of having to define four rates per additional site. While we still enforce the myosin head to interact with only one site at a time, it can now be any site with still only four rates to calibrate.

In (Chaintron, Caruel, and Kimmig 2023), we present this model in details using the concept of Poisson random measures3, and explain the conditions of its well posedness. We adapt the detailed balance condition (3.2) and show the ensuing compatibility with the second principle. We then propose a reduction strategy that allows to eliminate the variable \(h_t\). The resulting reduced model has only two degrees of freedom, namely \(s\) and \(\alpha\), like the Huxley-Hill model. The associated PDE formulation is analog to Equation 3.1 with a correction term. This shows that by extrapolation, more complex chemical-mechanical models (i.e. with more states) could be reformulated using this new framewotk. The paper also provide numerical illustrations of the h-model and of its reduced version.

3 The interested reader can find a short pedagogical introduction to the concept of Poisson Random measured in section 2.2 of (Chaintron et al. 2023).

3.4.2 Continuous model of the working stroke

Associated references

This wok was initiated by Marcucci and Truskinovsky (2010) and continued during M. Caruel’s PhD. Details about the model can be found in the following references (Caruel, Allain, and Truskinovsky 2013; Caruel and Truskinovsky 2018). Note also that (Caruel and Truskinovsky 2018), is already a review of several papers by Caruel and Truskinovsky.

The working stroke is the force-generating step of the actin-myosin interaction cycle, see Figure 1.3. It corresponds to a conformational change of the protein executed while it is attached to the actin filament. This conformational change can be interpreted as the swift relaxation towards thermal equilibrium of the internal degrees of freedom characterizing the conformation of the myosin protein, with an associated timescale of a few milliseconds (Huxley and Simmons 1971). This timescale is short compared to the time to complete the Lymnn-Taylor cycle (about 30 ms), which suggests that the power stroke is a purely mechanical phenomenon operating independently of ATP hydrolysis (Marcucci and Truskinovsky 2010; Caruel and Truskinovsky 2018).

In the classical framework of chemical-mechanical models presented in Section 3.2.1, the power stroke is viewed as a sequence of three of more chemical transitions (Huxley and Simmons 1971; Smith et al. 2008; Caremani et al. 2015; Månsson 2016). This representation is valid only in the limit where the energy wells characterizing these different states are “deep enough” for the Kramers escape rate approximation to be valid.

Figure 3.9: Continuous model of the working stroke. (a) Schematic representation of the cross-bridge model consisting of two degrees of freedom, the total elongation \(y\) and the variable \(x\) representing the conformation. (b) Conformational energy for chemical-mechanical framework of Huxley and Simmons (1971) (thin) or in the regularized continuous framework proposed by Marcucci and Truskinovsky (2010) (thick).

Based on the experimental observation, first made by Huxley and Simmons (1971), that the power stroke operates at fast timescales, Marcucci and Truskinovsky (2010) developed a model of the working stroke, where the conformational change is represented as a one-dimensional Langevin process in a double-well potential instead of a jump process, see Figure 3.9. This model can be considered as a regularized version of the original chemo-mechanical representation of the working stroke proposed by Huxley and Simmons (1971).

The model was introduced to model the collective response of groups of connected molecular motors to rapid load changes. This topic will be discussed in details in chapter Chapter 4. We now present how this continuous model of the power stroke can be coupled to a more classical “chemical”representation of the attachment-detachment process, to result in a hybrid jump-diffusion model actin-myosin interaction model.

3.4.3 Jump-diffusion model of the Lymn–Taylor cycle

Associated references

This work has been done in collaboration with the Inria MΞDISIM team. The formulation of the model and its first validation can be found in (Caruel, Moireau, and Chapelle 2019). The model was further refined and calibrated to reproduce cardiac muscle data during the PhD of F. Kimmig, see (Kimmig and Caruel 2020). Recently the formulation of the model was significantly improved, with the contribution of L.P. Chaintron, to streghthen its mathematical well-posedness and assess its compatibility with the thermodynamic principles (Chaintron et al. 2023).

Figure 3.10: Jump-diffusion model proposed by Chaintron et al. (2023). (a) The state of the myosin head is characterized by three stochastic variables: \(X_{t}\), \(Y_{t}\) and \(\alpha_{t}\), which represent the position of the head, the conformation and the state of attachment, respectively. The position of the nearest binding site is denoted by \(s\). The attachment-detachment process is characterized by four jumps whose rates are denoted \(k_{0\rightarrow 1}\), \(k^{\mathrm{rev}}_{0\rightarrow 1}\), \(k_{1\rightarrow 0}\) and \(k^{\mathrm{rev}}_{1\rightarrow 0}\). The energy landscape has two main ingredients: a quadratic term representing the elasticity of the head and a double-well potential associated with the conformational change. The latter potential depends on the attachment state. (b) illustration fo the Lymn-Taylor cycle with \(1\to2\): attachment; \(2\to 3\): power stroke; \(3\to 4\): detachment, and \(4\to 1\) recovery stroke. (a) is adapted from (Chaintron et al. 2023).

To extend the Marcucci and Truskinovsky (2010) model of the working stroke, and describe the whole Lymn–Taylor cycle, Caruel, Moireau, and Chapelle (2019) proposed a hybrid model, illustrated in Figure 3.10 (a). The state of the myosin motor is parametrized by the following stochastic variables

  • \(X_{t}\), the position of the myosin head with respect to its fixed location on the myosin filament,
  • \(Y_{t}\), the conformation of the head (power stroke),
  • \(\alpha_{t}\), the state of attachement (equal to 0 when the motor is detached and equal to 1 when the motor is attached).

As in the classical Huxley-Hill model, the myosin head can attach only to the nearest actin site, located at a distance \(s\) from the myosin anchor in the thick filament. When an attachment event (\(\alpha=0\to\alpha=1\)) occurs, the head position jumps to the value \(X_{t} = s(t)\).

The stochastic jump dynamics is formulated by Chaintron et al. (2023) using the concept of Poisson random measures, which allows to write the stochastic differential equation describing the dynamics of the discrete variable \(\alpha_{t}\). The jump process involves four transition rates

  • \(K_{0\rightarrow 1}(x,y,s)\) associated with the jump \(0 \rightarrow 1\) (direct attachment);
  • \(K^{\mathrm{rev}}_{0\rightarrow 1}(x,y,s)\) associated with the jump \(1 \rightarrow 0\) (reverse attachment);
  • \(K_{1\rightarrow 0}(x,y,s)\) associated with the jump \(1\rightarrow 0\) (direct detachment);
  • \(K^{\mathrm{rev}}_{1\rightarrow 0}(x,y,s)\) associated with the jump \(0\rightarrow 1\) (reverse detachment).

The conformation variable \(Y_{t}\) and the detached head position \(X_{t}\) both follow langevin dynamics. As in the classical Huxley-Hill formulation, the internal energy of the myosin motor depends on the attachment state. A typical actin-myosin interaction cycle is illustrated in Figure 3.10 (b).

Within this theoretical formalism, the complete jump-diffusion stochastic differential system can be written as

\[ \left\{ \begin{aligned} &\begin{split}\mathrm{d}\alpha_t &= \mathbb{1}_{\alpha_{t^-} = 0} \left[ \int_{\mathbb{R}_+} \mathbb{1}_{ z \leq K_{0\rightarrow 1}( X_{t^-} , Y_{t}, s(t) )} N_{0\rightarrow 1}( \mathrm{d}t , \mathrm{d}z ) \right.\\ &\left.\hspace{2cm}+\int_{\mathbb{R}_+} \mathbb{1}_{z \leq K^{\mathrm{rev}}_{1\rightarrow 0}( X_{t^-} , Y_{t}, s(t) ) } N^{\mathrm{rev}}_{1\rightarrow 0}( \mathrm{d}t , \mathrm{d}z ) \right]\\ & \phantom{ab}- \mathbb{1}_{\alpha_{t^-} = 1} \left[ \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} \mathbb{1}_{z \leq K_{1\rightarrow 0}( x, Y_{t}, s(t) )} N_{1\rightarrow 0}( \mathrm{d}t , \mathrm{d}x , \mathrm{d}z) \right.\\ & \hspace{2cm}+ \left.\int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} \mathbb{1}_{z \leq K^{\mathrm{rev}}_{0\rightarrow 1}( x, Y_{t}, s(t) )} N^{\mathrm{rev}}_{1\rightarrow 0}( \mathrm{d}t , \mathrm{d}x , \mathrm{d}z ) \right], \\ \end{split}\\ &\begin{split} \mathrm{d}X_t &= \mathbb{1}_{\alpha_{t^-} = 0} \left[ -\eta_x^{-1} \partial_x w_{0} ( X_t , Y_t ) \mathrm{d}t + \sqrt{2 \eta^{-1}_x k_B T}\, \mathrm{d}B^x_t \right] + \mathbb{1}_{\alpha_{t^-} = 1} \dot{x}_{c}(t) \mathrm{d}t \\ & \phantom{ab}+ \mathbb{1}_{\alpha_{t^-} = 0} \left[ \int_{\mathbb{R}_+} ( s({t}) - X_{t^-} ) \mathbb{1}_{z \leq K_{0\rightarrow 1}( X_{t^-} , Y_{t}, s(t) )} N_{0\rightarrow 1}( \mathrm{d}t , \mathrm{d}z ) \right.\\ & \hspace{2cm}+ \left.\int_{\mathbb{R}_+} ( s(t) - X_{t^-} ) \mathbb{1}_{z \leq K^{\mathrm{rev}}_{1\rightarrow 0}( X_{t^-} , Y_{t}, s(t) ) } N^{\mathrm{rev}}_{1\rightarrow 0}( \mathrm{d}t , \mathrm{d}z ) \right] \\ & \phantom{ab}+ \mathbb{1}_{\alpha_{t^-} = 1} \left[ \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} ( x - s(t) ) \mathbb{1}_{ z \leq K_{1\rightarrow 0}(x, Y_{t}, s(t) )} N_{1\rightarrow 0}( \mathrm{d}t , \mathrm{d}x , \mathrm{d}z) \right. \\ & \hspace{2cm}+ \left.\int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} ( x - X_{t^-} ) \mathbb{1}_{z \leq K^{\mathrm{rev}}_{0\rightarrow 1}(x, Y_{t}, s(t) )} N^{\mathrm{rev}}_{0\rightarrow 1}( \mathrm{d}t , \mathrm{d}x , \mathrm{d}z) \right], \\ \end{split}\\ &\mathrm{d}Y_t = -\eta_y^{-1} \partial_y w_{\alpha_t} ( X_t , Y_t ) \mathrm{d}t + \sqrt{2 \eta^{-1}_y k_B T}\, \mathrm{d}B^y_t, % \end{split} \end{aligned} \right. \tag{3.4}\]

where \(\eta_{x,y}\) are drag coefficients, and \((B_{t}^{x})_{t\geq 0}\) and \((B_{t}^{y})_{t\geq 0}\) are independent Brownian motions. The meaning of the integrals of the form \(\int_{0}^{t}\int_{\mathbb{R}_{+}}F(t,z)N(\mathrm{d}t,\mathrm{d}z)\) is explained in more details with extra references and numerical implementation in (Section 2.2 and 4.1 of Chaintron et al. 2023). Notice that, in the second equation, the motion of \(X_{t}\) is determinitic when the head is attached, like in the two state model presented in Section 3.2.1.

The peculiarity of this model lays in the assymetry between the attachment and the detachment jumps, which results from the different system’s dimensionality in the attached and detached states. When the motor is detached, it evolves in a two dimenional space \((X_{t},Y_{t})\), while when it is attached, \(X_{t}=s(t)\) so the dynamics besomes one dimensional. Hence, when an attachment jump occurs from a random configuration \((X_{t},Y_{t})\), the dynamics is then projected on a one dimensional space, parametrized by \(s(t)\), where the value of \(X_t\) after the jump is fully determined by the position of the site.

The situation is different with the detachment, since there is no a priori value prescribed for \(X_{t}\) after the jump. This discrete-to-continuous detachment jump must then select a radom location \(x\) for the detached myosin head. The probability law that choose this new location have to be defined carefully to preserve the thermodynamic compatibility.

This problem was solved in (Chaintron et al. 2023) by considering the total jumps \((0,x,y,s)\to (1,s,y,s)\to(0,x,y,s)\) and writing the detailed balance conditions \[ \left\{ \begin{aligned} K_{0\rightarrow 1}(x,y,s) &= h_x \exp \biggl[ -\frac{w_1(s,y) - w_0(x,y) }{k_B T} \biggr] K^{\mathrm{rev}}_{0\rightarrow 1}(x,y,s),\\ K_{1\rightarrow 0}(x,y,s) &= h^{-1}_x \exp \biggl[ -\frac{(w_0(x,y) - \mu_T) - w_1(s,y) }{k_B T} \biggr] K^{\mathrm{rev}}_{1\rightarrow 0}(x,y,s), \end{aligned} \right. \tag{3.5}\] where \(h_{x}\) is a characteristic length for \(x\), taken equal to the characteristic length of the power stroke. Provided Equation 3.5, only two rates have to be calibrated. From the rate \(K_{1\rightarrow 0}\) one can define the global detachment rate, describing the frequency of the detachment, considering all detached positions \(x\) that can be reached from an attached configuration \((1,s,y,s)\) \[ k_{1\rightarrow 0}(y,s) = \int_{-d/2}^{d/2}K_{1\rightarrow 0}(x,y,s)\mathrm{d}x, \] A similar rate \(k^{\mathrm{rev}}_{0\rightarrow 1}\) can be defined for the reverse attachment jump rate \(K^{\mathrm{rev}}_{0\rightarrow 1}\). The rates \(k_{1\rightarrow 0}\) and \(k^{\mathrm{rev}}_{0\rightarrow 1}\) are analog to the rates presented in Section 3.2.1. In the numerical implementation of the process, the global rates \(k_{1\rightarrow 0}\) and \(k^{\mathrm{rev}}_{0\rightarrow 1}\) are used to termine the time of the jump. When a detachment jump occurs, the jump amplitude is drawn from the probability law \(k_{1\rightarrow 0}(y,s)^{-1}K_{1\rightarrow 0}(x,y,s)\mathrm{d}x\) (or \(k^{\mathrm{rev}}_{0\rightarrow 1}(y,s)^{-1}K^{\mathrm{rev}}_{0\rightarrow 1}(x,y,s)\mathrm{d}x\)), see (Section 4.1 and 4.2 of Chaintron et al. 2023) for the details.

Figure 3.11: Hill force-velocity curve produced by the stochastic jump-diffusion model of Chaintron et al. (2023). The curve is obtained by computing the average force of 105 realizations in the steady state regime under imposed sliding velocity. (a) steady state force \(T_{c}\) (normalized to isometric force \(T_{0}\)) as a function of the shortening velocity \(-\dot{x}_{c}\). (b) distribution of the system in the attached state. The level sets reproduce the energy landscape. the colored dots show three snapshots of the stochastic realizations corresponding to the points A, B and C shown in (a) for a subset of \(5\times 10^{3}\) randomly chosen realizations. (c) Distributions of attached heads at points A, B and C. (d) Distribution of the conformational variable \(Y_{t}\) in the double well potential for selected relative position \(s\) [marked * beside the vertical lines in (b)] of the binding site.

The model developed in Chaintron et al. (2023) has been calibrated to reproduce data obtained on rat trabeculae in various experimental studies. For the details about the calibration procedure, we refer to (Kimmig and Caruel 2020; Chaintron et al. 2023). In Figure 3.11, we show an extract of the results illustrating the Hill force-velocity curve. The model faithfully reproduces the fast transient response to rapid load changes (see figure 6 and 7 of Chaintron et al. 2023), and can be used to predict the thermodynamic efficiency of the contraction.

3.5 Ongoing work

3.5.1 Molecular dynamics study of the Lymn–Taylor cycle

Note

This topic is the motivation of R. Manevy’s the PhD thesis at MSME. It is done in collaboration with A. Houdusse from Institut Curie. This work has been presented in several conferences and workshop (Caruel et al. 2022; Manevy, Detrez, et al. 2021; Manevy, Caruel, et al. 2021). A paper is in preparation.

The objective of this work is to close the gap between structural mechanisms of the molecular motors functioning inferred from crystallographic data and more coarse theoretical models. The strategy is to project the high-dimension dynamics of the molecular system onto a low dimension manifold, enabling physiological timescale simulations based on Stochastic Differential Equations alike the ones presented in section Section 3.4.3.

For instance, the parameters characterizing energy landscape that is used in (Marcucci and Truskinovsky 2010; Caruel, Moireau, and Chapelle 2019; Kimmig and Caruel 2020) to describe the power stroke are usually calibrated using data from fast load change experiments performed on muscle fiber. One could also compute this energy landscape from molecular dynamics simulations.

In the PhD work of R. Manevy, we studied the step of the actin-myosin interaction cycle where the inorganic phosphate (Pi), a product of the ATP hydrolysis, is released from the active site. The interplay between this escape from the active site and the power stroke conformational change is a subject of debate: is the release occuring before or after the working stroke, or can they occur simultaneously (see Llinas et al. 2015 for a review on the topic)?

To test several escape pathways and their associated energetic costs, we performed Umbrella Sampling simulations starting with two christallographic structures and several local configuraiton of the Pi within the active site, see Figure 3.12.

Figure 3.12: Active site of the Myosin VI structures pre-powerstroke state (PPS) (PDB 2V26) (A and B) and Pi-release state (PiR) (PDB 4PJM) (C and D). The P-Loop (residues 150 to 157) in brown, Switch I (residues 195 to 205) in pink and Switch II (residues 456 to 467) in orange are represented in cartoon. Other subdomains of the structure are not represented (see Figure 3.3). Pi, ADP, Mg2+, ARG199 and ARG205 of Switch I, GLU461 of Switch II are represented in licorice with color code C cyan, H white, O red, N blue, P brown, and Mg2+ light pink. In each panel, Switch II is represented both in solid transparent to show its displacement between the PPS and the PiR conformation. The Pi release pathways observed during our US simulations are represented by red double arrows following the terminology of Cecchini, Alexeev, and Karplus (2010). The Switch II displacement is characteristic of the opening and closing of the back door I escape pathway for Pi.

The constraint is the distance between the Pi and the magnesium ion associated with ADP in the active site. By increasing this distance in a step-by-step manner, we force the inorganic phosphate out of the active site. The potential of mean force associated with this process can be reconstructed, giving an estimate of the energy barrier associaed with the escape pathway taken by the Pi. An example of the escape process is illustrated in Figure 3.13. Only the strutures surrounding the active site are shown for clarity.

Figure 3.13: Umbrella sampling simulation of the escape of Pi from the active site. The video shows the interactions between the nucleotides binding loops during the umbrella sampling simulation for different protonation states of the Pi (named PiR A to E). The simulations are displayed from the point of view of the back door I. The color code is similar to Figure 3.12. Courtesy of F. Detrez (MSME). Data produced by R. Manevy.

We found that the initial condition of the simulation, and in particular the orientation of the Pi within the active site, greatly affects the escape pathway and the associated potential of mean force. The statistical validation of these observations are work in progress.

3.5.2 Purely mechanical model of the actin-myosin interaction

Note

This work is a follow-up of our collaboration with L. Truskinovsky and R. Sheshka. Preliminary work can be found in (Sheshka and Truskinovsky 2014) and (section 4.2 of Caruel and Truskinovsky 2018).

In Section 3.4.3, we presented a hybrid chemical-mechanical framework of the actin-myosin interaction. As mentioned in Section 3.2, there exists a category of models where the dynamics of the molecular motor is represented as a Langevin-type drift-diffusion process with a colored noise. The prototypical example of this class of models was proposed by Magnasco (1993), see Equation 3.3.

Figure 3.14: Illustration of the stochastic mechanical model proposed by Sheshka and Truskinovsky (2014), see also Caruel and Truskinovsky (2018). The interaction between actin and myosin is represented by a periodic potential \(\Phi\). The position of the lever arm of the protein (\(y(t)\)) is sumbitted to a periodic force \(f(t)\) that mimicks the energy input from ATP hydrolysis.

The principle of this type of model is illustrated in Figure 3.14. The molecular motor is parametrized by the position of the head \(X_{t}\) and the position of the tip of the lever arm \(Y_{t}\). Both variables are endowed with Langevin stochastic dynamics. The associated potential includes (see Figure 3.15)

  • the energy landscape \(\Phi\) which represents the interaction between the tip of the head and the actin track (a),
  • the bistable potential \(u_{\text{ss}}\) representing the power stroke conformational change energetics (b).

Figure 3.15: Ingredients of the stochastic mechanical model proposed by Sheshka and Truskinovsky (2014). (a) actin interaction potential \(\phi\); (b) power stroke potential \(u_{\text{ss}}\); (c) colored component of the noise taking the form of a periodic force acting on the variable \(y-x\); (d) representation of biochemical-mechanical coupling through the variable \(d\) representing the affinity of the myosin head to actin, depending of the power stroke state.

As in the Magnasco (1993) model, the colored component of the noise is a piecewise constant periodic force \(f\) that operates on the variable \(Y_{t}-X_{t}\), i.e. on the power stroke, see Figure 3.14 (c). The action of this force on the power stroke is the mechanical representation of the energy input provided by ATP hydrolysis. The directionality of the motion can be achieved either by making the actin interaction potential \(\Phi\) non symetric or by introducing a control mechanism \(d(Y_{t}-X_{t})\in[0,1]\), that mimicks the coupling between mechanics and biochemistry in the Lymn–Taylor cycle. The potential of the system can then be written as \[ G(x,y,d,t) = d\Phi(x) + u_{\text{ss}}(y-x) + f(t)(y-x) - f_{\text{ext}}y, \tag{3.6}\] where \(f_{\text{ext}}\) is the external force. The dynamics of the system is then written as \[ \left\{ \begin{aligned} \mathrm{d}X_{t} &= -\partial_{x}G(X_{t},Y_{t},d,t)\mathrm{d}t + \sqrt{2D}\mathrm{d}B_{t}^{x},\\ \mathrm{d}Y_{t} &= -\partial_{y}G(X_{t},Y_{t},d,t)\mathrm{d}t + \sqrt{2D}\mathrm{d}B_{t}^{y}. \end{aligned} \right. \] In Equation 3.6, the coupling \(d\) characterizes the strength of the attachment. When \(d=1\), the actin potential affects on the evolution of \(X_{t}\), effectively modeling the attachment. When \(d=0\), the interaction with actin is turned “off”, which represents the detached state.

Considering the dependence of \(d\) on the power stroke conformation \(Y_{t} - X_{t}\) shown in Figure 3.15 (d), the multiplicative combination of \(d\) and \(\Phi\) effectivety mimicks the loss of affinity of the myosin head for actin as the power stroke conformational change happens. At the begining of the power stroke \(Y_{t}-X_{t}\approx -a/2\), so \(d=1\). Therefore, the gradient actin potential participates in driving the position \(X_{t}\). After the completion of the power stroke \(Y_{t}-X_{t}\approx a/2\), so \(d=0\). Therefore, the contribution of the actin potential to the motion of the head vanishes, mimicking a detached state.

An interesting characteristic of this family of model involving a control parameter is that directionality can be achieved even if all the potentials are symetric, like shown in Figure 3.15. More refined version of the coupling \(d\) can incorporate memory effect which can, according to the preliminary results of Sheshka and Truskinovsky (2014), augment the mechanical performance of the motor, see Figure 3.16.

Figure 3.16: Particle trajectory produced in a load-free setting by a mechanical stochastic model with biochemical-mechanical coupling including memory effect. The coupling parameter \(d\) dependence on the power stroke state \(y-x\) takes the form of a hysteresis (see right column). Changing the width of this hysteresis has a strong impact on the steady state velocity of the molecular motor. Figure taken from (Sheshka and Truskinovsky 2014).

3.6 Other project: nano diffusion in the bone tissue

Associated references

The results of this work are published in (Prakash, Lemaire, Di Tommaso, et al. 2017; Prakash, Lemaire, Caruel, et al. 2017).

The biomechanics team of the MSME laboratory actively work on the mechanics of bone remodeling at various scales. Biochemical signaling for bone remodeling involves ionic diffusion within nanopores that connect bone cells to each other. We used molecular dynamics tools to simulate confined diffusion at the nanoscale and determine the physical parameters (diffusion tensor, viscosities, etc.) associated with the flow in the nanopores. This work enabled the determination of relevant parameters for a multiscale model of bone remodeling and also contributed to the understanding of the molecular mechanism of intercellular signaling within the bone tissue. For more details we refer to (Prakash, Lemaire, Di Tommaso, et al. 2017; Prakash, Lemaire, Caruel, et al. 2017).

3.7 Remaining challenges and future work

  • A natural follow-up of the work on the jump-diffusion model presented in Section 3.4.3 is to reconsider the mean-field hypothesis and model a finite size cluster of interacting molecular motors. This line of research was started in Caruel, Allain, and Truskinovsky (2013) and will be presented in more detailed in the next chapter.
  1. The question raised is the representation of the mode of interaction between the motors. One approach consists in considering short range elastic interactions between motors that are attached on the same actin filament. In this way, next-to-nearest neighbors cooperative mechanisms can emerge and propagate inside the contractile unit.
  2. This extension of a finite size system is useful to compare the model predictions with experimental measurements performed on synthetic nanomachines mimicking the muscle contractile units, see (Pertici et al. 2018).
  • The potential next steps of R. Manevy’s thesis discussed in Section 3.5.1 are to extend the current study to the other steps of the Lymn–Taylor cycle and to strengthen of the methodology used for determining the adequate collective variables characterizing the conformational changes.

  1. This work heavily relies on the availability of crystallographic structures serving as initial (or final) condition for these simulations. These structures are not always available in the protein databank, which limits the applicability of a model calibration approach based solely on molecular dynamcis simulations. In particular, it still seems difficult to simulate detachment or attachments events using these techniques before establishing the structural basis of the affinity of myosin for actin.

  2. Choosing appropriate collective variables to describe the motion of the protein during a conformational change is not easy. The energy landscape characterizing the protein is high-dimensional and very rugged so that (i) the collective variables found numerically may not be the most representative (i.e. the fastest) and (ii) the sampling of the effective free energy landscape ineffective. Our future work will involve exploration of various sampling methods and potential use of AI based techniques, see (Zhang et al. 2023).

Back to top