First-passage times in complex energy landscapes: a case study with nonmuscle myosin II assembly

In order to give an instructive example of a complex one-dimensional energy landscape with many potential barriers of different height, we first introduce the case of the self-assembly of NM2-molecules into minifilaments. As schematically depicted in figure 1(a), one NM2 monomer consists mainly of two identical myosin II heavy chains (shown in light and dark green) that dimerize to form its rigid helical backbone. Hereby, each heavy chain begins with a globular motor head domain at the N-terminal containing bindings sites for actin filaments as well as for ATP-molecules. Following the head region comes a neck region in which two regulatory light chains and two essential light chains bind (shown in blue). Then comes the roughly 160 nm long and relatively rigid coiled-coil rod and finally the non-helical tailpiece that disrupts the coiled-coil motif. Although the whole complex is in fact a hexamer, here we call it a NM2-molecule for simplicity. In the coiled-coil rod domain the two heavy chains wind around each other with an axial spacing of 0.1456 nm between adjacent residues and a pitch of 3.5 residues per turn due to interactions between periodically placed hydrophobic residues [102, 103]. For this purpose, each NM2 heavy chain features an amino acid sequence with a heptad repeat labeled from a to g as schematically depicted in the inset of figure 1(a). Here, residues at positions a and d form the hydrophobic interface between two adjacent side chains, whereas the remaining residues are hydrophilic and exposed to the solvent. The assembly of the NM2-molecules into dimers is mediated by the electrostatic interactions of the exposed residues, whose charge pattern leads to either parallel or antiparallel staggers, as shown in figure 1(b). The predominance of electrostatic driving forces during the assembly process was demonstrated by in vitro experiments where an increase in salt concentration was found to inhibit filament formation due to electrostatic screening effects [104]. The final minifilament contains up to 30 NM2-molecules, is bipolar and around 300 nm large, compare figure 1(c).

Zoom In Zoom Out Reset image size

The derivation of the electrostatic interaction energy between two NM2-molecules relies on the procedures recently introduced in reference [97]. First the amino acid sequence of the coiled-coil rod domain of individual NM2-molecules, as provided by the NCBI protein database [105], is translated into a linear chain of point-charges. We have analyzed the amino acid sequences of the three human NM2-isoforms (denoted as NM2A, NM2B and NM2C) via the sequence-based coiled-coil prediction software Paircoil2 [106]. By employing pairwise residue probabilities [107], Paircoil2 is able to identify the location of coiled-coil motifs in protein amino acid sequence data. Performing this statistical analysis yields a p-value for each individual residue of the sequence data, which in turn allows us to test the null hypothesis whether the respective residue is not located within a coiled-coil domain. As depicted in figure 2(a), this hypothesis is rejected for a domain ranging from roughly the 850th up to the 1920th amino acid residue for all NM2-variants, indicating that these residues are indeed part of the coiled-coil motif of the NM2 heavy chain. The first large non coiled-coil region indicates the motor and neck region, whereas the small region at the end of the amino sequence represents the non-helical tailpiece (compare figure 1(a)).

Zoom In Zoom Out Reset image size

After identifying the NM2 coiled-coil motif ($\sim 1100$ out of total $\sim 2000$ amino acids) and the positions of individual residues along the rod, we are now able to associate all required amino acids of the coiled-coil with their corresponding charges. Amino acids at residue position a and d are hereby assigned zero net charge as they face inwards and form the hydrophobic coiled-coil backbone (compare inset in figure 1(a)). The remaining residue positions are either assigned a charge of '+2e' for positively charged amino acids (i.e. arginine and lysine) or '−2e' for negatively charged amino acids (i.e. glutamic acid and aspartic acid). Moreover point-like charges are arranged with an axial spacing of Δx = 0.1456 nm between adjacent neighbors [97, 102, 103] which subsequently yields the correct dimensions for the NM2-molecule length of approximately 160 nm. The corresponding distribution of charges along the NM2 coiled-coil rod domain is depicted in figure 2(b) for all three different NM2-variants. For clarity, we hereby smooth the distribution of charges via a sliding window of 98 residues [97, 99]—the preferred axial shift between interacting molecules due to the 196 residue repeat—which allows us to identify the positive C-terminal tip assembly critical domain (ACD), and the five regions with negative charge spaced along the rod domain. We note that the quantitative difference between the different NM2-variants is comparatively small.

The + − − − − − charge pattern revealed in figure 2(b) is the physical basis of the self-assembly of minifilaments and has been identified before by similar procedures [97–99]. In line with these previous modeling approaches, our calculations of energy landscapes are restricted to the electrostatic interactions of two aligned NM2-molecules which are based on the charge distribution calculated in section 2.1. For simplicity we do not account for bending of NM2 rods that has been used before by including bending energies to complement the electrostatic interactions [97] as it would require additional numerical procedures and effectively renders the energy landscapes two-dimensional. We also neglect potential contributions of the entropy of segment bending to the free energy. Because the persistence length of the myosin II rods is 130 nm [108, 109] and dangling parts to the polymers are usually in the same order of magnitude or less, such contributions are expected to contribute only little in our context. To calculate the electrostatic interaction energy between two rods, we note that under physiological salt concentrations we are in the weak coupling regime for a rod with 0.085 e per residue, a distance 0.1456 nm between residues and a radius 1 nm [97]. Thus the appropriate treatment of electrostatics is Poisson–Boltzmann theory [110, 111]. Our considered system is additionally characterized by small ionic strengths, allowing us to linearize the Poisson–Boltzmann equation, yielding the Debye–Hückel equation

$$\begin{equation}{\nabla }^{2}\psi (\mathbf{r})={\kappa }^{2}\psi (\mathbf{r}).\end{equation} \tag{ 1 }$$

Here ψ denotes the electrostatic potential and κ is the inverse Debye–Hückel screening length that sets the typical length scale of the electrostatic screening.

To now compute the electrostatic interaction energy of two aligned NM2 rods, each modeled as a linear chain of charges, we always have to consider two different cases. As schematically illustrated in figure 1(b), NM2-molecules assemble in either parallel or antiparallel fashion, and they may align with a different stagger s. For the geometry of a single point-particle with charge q, it is possible to analytically solve equation (1), resulting in a Yukawa-type potential [110] of the form

$$\begin{equation}\psi (r)=\frac{q}{4\pi {\epsilon}{{\epsilon}}_{0}}\frac{\mathrm{exp}(-\kappa r)}{r}.\end{equation} \tag{ 2 }$$

Exploiting the linearity of the Debye–Hückel equation (1) furthermore allows us to directly compute the total electrostatic interaction energy U(s) of two aligned NM2-rods as a function of their stagger s as

$$\begin{equation}U(s)=\sum\limits _{i=1}^{N}\sum\limits _{j=1}^{M}\frac{{q}_{i}{q}_{j}}{4\pi {\epsilon}{{\epsilon}}_{0}}\frac{\mathrm{exp}(-\kappa {r}_{ij})}{{r}_{ij}},\end{equation} \tag{ 3 }$$

where we define ${r}_{ij}\equiv \sqrt{{y}^{2}+{[(i-j){\Delta}x-s]}^{2}}$ as the distance between residues. Point charges are denoted with q whereas N and M denote the total number of residues of the respective rods with i and j as their corresponding index. The distance along the NM2 rod, i.e. the main axis, is given by x, while y denotes the lateral spacing between two aligned NM2-molecules. In the case of an antiparallel NM2-dimer configuration we invert one of the two charge sequences.

The resulting electrostatic interaction energy landscapes for NM2-dimers are shown in figures 3(a) and (b) as a function of the stagger s and in units of the thermal energy k_B T = 1/β for parallel and antiparallel configurations, respectively. By calculating the interaction energy according to equation (3) we assume a lateral spacing of y = 2 nm [87, 97]. Moreover, we do not take negative staggers into account as their positive ACDs would not be in contact with the negative charged region, i.e. their interaction would become highly unfavorable. The experimentally known staggers (s = 14.3 nm, s = 43.2 nm and s = 72.0 nm for parallel alignment and s = 113 nm to s = 118 nm for antiparallel alignment) are marked in figures 3(a) and (b) by black arrows and emerge as local minima generated by the electrostatic interaction between the positive ACD and regions with increased net negativity along the NM2-rod. However, as noted earlier, these minima are not very clear, which motivates our approach to use these energy landscapes to investigate the corresponding FPT-problems.

Zoom In Zoom Out Reset image size

To investigate situations in which the NM2-dimer is being pulled apart, e.g. due to contractile forces in the actomyosin network, we study the impact of an external pulling force F on the energy landscape. Force-mediated rearrangements may be of particular importance at the beginning stages of minifilament assembly and during the dynamical partitioning of minifilaments [91, 112], when filaments physically split into a pair of two 'daughter' filaments. To account for the mechanical load we apply an additional constant external pulling force F > 0 pN along the stagger main-axis coordinate in positive s-direction. As seen in figures 3(c) and (d) for parallel and antiparallel alignment, respectively, this driving leads to a tilted NM2-dimer energy landscape

$$\begin{equation}\beta {U}^{(F)}(s)=\beta {U}_{0}(s)-\beta Fs,\end{equation} \tag{ 4 }$$

where U₀ ≡ U^(0 pN)(s) denotes the electrostatic potential energy as given by equation (3) of the original 'untilted' NM2-dimer.

Our numerical treatment of equation (10), which is based on continuous-state dynamics of an overdamped Fokker–Planck equation, has shown a pattern in the MFPT as a function of stagger that requires an explanation in terms of the key features of the underlying complex energy landscape. To identify these features, we now coarse-grain the continuous-state dynamics into a set of discrete states to further analyze the MFPT as a function of the initial stagger and learn more about the connection between energy landscape and the shape of the MFPT. The motivation hereby is twofold. First, the existence of high energy barriers that are significant to the thermal energy k_B T (see figure 3) allows us to motivate a clear time-scale separation such that a diffusing particle would spend the majority of its time located in valleys of the energy landscape before eventually escaping the barrier—a situation well described by stochastic transitions on a discrete network in terms of a Markovian jump process. Second, the MFPT until a NM2-dimer detaches has, compared to the underlying energy landscape, a relatively simple shape suggesting that only a selected few potential barriers and minima contribute to the overall average time needed. As seen in figure 4 the shape of the MFPT is characterized by a of number pronounced decreases which coincide with the position of specific potential barriers. Unfortunately this feature remains elusive for the continuous-state description as the expression used to calculate the MFTP ${\mathcal{T}}_{1}$ in equation (10) poses a non-trivial transformation of the underlying energy landscape.

Within the coarse-graining procedure we treat each local minima of the corresponding NM2 energy landscapes as a discrete state in the network description. Being in a state is therefore equivalent to residing in the associated local basin. By splitting the energy landscape into a number of N basins (compare figure 6(b)), delimited by their adjacent potential maxima, the resulting Markovian network is then described by the set Ω = {0, 1, ..., N − 1} of N discrete states. Local minima in the energy landscape are identified numerically using the signal.argrelmin function of the Python library SciPy [116]. The network dynamics follow the master equation

$$\begin{equation}{\partial }_{t}{p}_{i}(t)=\sum\limits _{j=1}^{N-1}[{w}_{ji}{p}_{j}(t)-{w}_{ij}{p}_{i}(t)].\end{equation} \tag{ 11 }$$

Here p_i (t) represents the occupation probability to be in state i at time t and w_ij denotes the transition rate to pass from state i to state j. Transitions only occur between directly neighboring states due to the considered problem of stagging two neighboring rods. Therefore, the master equation (11) may be re-written in a one-step form as

$$\begin{equation}{\mathrm{\partial }}_{t}{p}_{i}(t)={\lambda }_{i-1}{p}_{i-1}(t)+{\mu }_{i+1}{p}_{i+1}(t)-({\lambda }_{i}+{\mu }_{i}){p}_{i}(t),\end{equation} \tag{ 12 }$$

where we define the transition rates λ_i = w_i,i+1 and μ_i = w_i,i−1 (transition rates between non-neighboring states vanish). Figure 6(a) shows the corresponding reaction network. The MFPT is now identified with the time at which the transition N − 2 → N − 1 into the absorbing state takes place. As before, the state N = 0 is reflecting.

Zoom In Zoom Out Reset image size

Without loss of generality, we now choose the NM2A dimer in antiparallel configuration as an exemplary case. Figure 6(b) shows the N = 39 states that we obtain numerically from the corresponding energy landscape. For Markovian network dynamics, the exit rate r_i (or inverse lifetime) of the state i, i.e., the time that passes until the state changes, is exponentially distributed according to ${p}_{i}^{\text{exit}}(t)={r}_{i}\,\mathrm{exp}(-{r}_{i}t)$ with ${r}_{i}=1/{\langle t\rangle }_{i}^{\text{exit}}={\sum }_{j\in {\mathcal{N}}_{i}}{w}_{ij}$ (here ${\mathcal{N}}_{i}$ labels the neighbors of state i). The exit time distribution is independent of the final state j, allowing us to write the (splitting) probability of the transition from state i to j (i.e., that the next visited state will be j conditioned on that the system is in state i) as ϕ_j|i = w_ij /r_i . Under these considerations the transition rates are then obtained as

$$\begin{equation}{w}_{ij}=\frac{{\phi }_{j\vert i}}{{\langle t\rangle }_{i}^{\text{exit}}}.\end{equation} \tag{ 13 }$$

Quantities like the average exit time and the splitting probability entering equation (13) can be related to a conditional FPT problem [78, 117]. Here, one seeks the first time a system reaches a specific target state conditioned on that it has not yet reached any other (of possibly arbitrary many) target states before. Based on this framework and using the renewal theorem [118], results for the splitting probability and the mean exit time have been derived recently for general networks with a 'star-like' topology [117]. The 'inner' node is thereby the starting state and all 'outer' nodes are target states. A line network as we consider it here is a particular simple case of this network topology. Hereby, assuming that the full network dynamics between two network states evolve according to a continuous space-time overdamped Markovian diffusion, the splitting probability ϕ_j|i that starting from some state i the neighboring state visited next will be state j is given by the explicit expression

$$\begin{equation}{\phi }_{j\vert i}=\frac{1/{I}_{j\vert i}^{(1)}}{{\sum }_{k\in {\mathcal{N}}_{i}}1/{I}_{k\vert i}^{(1)}}\end{equation} \tag{ 14 }$$

and is normalized such that ${\sum }_{j\in {\mathcal{N}}_{i}}{\phi }_{j\vert i}=1$. The second important quantity required to obtain transition rates is the mean exit time of state i which is obtained using the same framework via the expression

$$\begin{equation}{\langle t\rangle }_{i}^{\text{exit}}=\sum\limits _{k}\;{\phi }_{k\vert i}{I}_{k\vert i}^{(2)},\end{equation} \tag{ 15 }$$

where the two auxiliary integrals ${I}_{j\vert i}^{(1)}$ and ${I}_{j\vert i}^{(2)}$ in equations (14) and (15) are defined [117] as

$$\begin{equation}{I}_{j\vert i}^{(1)}={\int }_{0}^{{l}_{j\vert i}}\mathrm{d}{y}_{1}\,{\ g}_{j\vert i}^{(1)}\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad {g}_{j\vert i}^{(1)}=\frac{1}{{D}_{j\vert i}}\,{\text{e}}^{\beta {U}_{j\vert i}({y}_{1})},\end{equation} \tag{ 16 }$$

$$\begin{equation}{I}_{j\vert i}^{(2)}={\int }_{0}^{{l}_{j\vert i}}\mathrm{d}{y}_{1}{\int }_{0}^{{y}_{2}}\mathrm{d}{y}_{2}\,{\ g}_{j\vert i}^{(2)}\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad {g}_{j\vert i}^{(2)}=\frac{1}{{D}_{j\vert i}}\,{\text{e}}^{\beta {U}_{j\vert i}({y}_{1})-\beta {U}_{j\vert i}({y}_{2})}.\end{equation} \tag{ 17 }$$

Here l_j|i denotes the distance between two adjacent minima (i.e. states i and j) defined through the underlying energy landscape and U_j|i is the local energy landscape relative to the value at state i between states i and j as schematically depicted in figure 6. For the purpose of calculating the transitions rate between two neighboring states we therefore identify their intervening local potentials and numerically compute the integrals in equations (16) and (17) employing the trapezoid rule function integrate.trapezoid as implemented in the python library SciPy [116].

After identifying all network states and calculating the corresponding transition rates w_ij , and thus λ_i and μ_i of equation (12), the Markov jump process is now fully specified. Therefore, we are able to determine the MFPT, i.e. the average detachment time of NM2-dimers, within the coarse-grained network description of the minima-to-minima dynamics of a diffusing Brownian particle in the original overdamped Fokker–Planck picture. As for the continuous system, the corresponding absorbing state is hereby always placed in the last state such that reaching it corresponds to the detachment of the NM2-molecules.

In analogy to the analytical solution of the MFPT in equation (10), which holds for a Fokker–Planck equation for overdamped Langevin dynamics in an external potential U, for a Markov process confined between mixed boundaries (a reflecting and b absorbing with a < b) the MFPT is given by [78, 119]

$$\begin{equation}{\mathcal{T}}_{1}({i}_{0};a,b)=\sum\limits _{i={i}_{0}}^{b-1}\sum\limits _{n=a}^{i}{\lambda }_{n}^{-1}{{\Pi}}_{n+1,i}\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad {{\Pi}}_{i,j}=\frac{{\mu }_{i}{\mu }_{i+1}\dots {\mu }_{j}}{{\lambda }_{i}{\lambda }_{i+1}\dots {\lambda }_{j}}\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace i\leqslant j,\end{equation} \tag{ 18 }$$

where i₀ is the starting state, a the reflecting state, and b the absorbing target state. In our case this MFPT has to be calculated as a function of the initial stagger i₀ and a and b are set to the first and last state, respectively. This calculation now is much more computationally efficient than the continuum version. Figures 7(a) and (b) show a direct comparison between the MFPT ${\mathcal{T}}_{1}$ computed using equation (10) (lines) and equation (18) (symbols) for the NM2A-dimer. Clearly, the results of the two different approaches agree very well with each other, thus validating our coarse-graining approach. The same picture emerges for the other NM2 variants. We additionally plot the mean exit-time ${\langle t\rangle }_{i}^{\text{exit}}$ from equation (15) as bars. Interestingly, the experimentally observed staggers (black arrows) for the parallel configuration corresponds to large values for ${\langle t\rangle }_{i}^{\text{exit}}$. In particular, the most prominent stagger at s = 14.3 nm corresponds to the highest overall mean exit-time. While the mean exit-times corresponding to experimentally observed staggers for antiparallel alignment are not notably large on a global scale, they are indeed large compared to adjacent states. This suggests that these staggers confer some local robustness against perturbations. We additionally note that the periodicity between small mean exit-times seems to follow the characteristic 98-residue repeat pattern already encountered earlier.

Zoom In Zoom Out Reset image size

The relatively simple shape of the MFPT and the characteristic feature to decrease substantially over several orders of magnitude by only a slight variation of the initial stagger position between two NM2-molecules in the detachment process remained elusive in the continuous description. The discrete nature of the Markov jump process now allows us to gain further insights in the behavior of the MFPT of the detachment process of two aligned NM2-molecules. We first note that an exploration of the complete network topology (i.e. the complete energy landscape) takes place if the system is initially prepared in state i₀ = 0. Then is has to visit each state at least once before reaching the absorbing target at the other end of the network chain. For this case equation (18) is explicitly given by

$$\begin{align}\hfill {\mathcal{T}}_{1}(0;0,b)& =\sum\limits _{i=0}^{b-1}\sum\limits _{n=0}^{i}{\lambda }_{n}^{-1}{{\Pi}}_{n+1,i}=\underset{i=0}{\underbrace{{\lambda }_{0}^{-1}}}+\underset{i=1}{\underbrace{{\lambda }_{0}^{-1}{{\Pi}}_{1,1}+{\lambda }_{1}^{-1}}}+\underset{i=2}{\underbrace{{\lambda }_{0}^{-1}{{\Pi}}_{1,2}+{\lambda }_{1}^{-1}{{\Pi}}_{2,2}+{\lambda }_{2}^{-1}}}\hfill \\ \hfill & \quad +\cdots +\underset{i=b-1}{\underbrace{{\lambda }_{0}^{-1}{{\Pi}}_{1,b-1}+{\lambda }_{1}^{-1}{{\Pi}}_{2,b-1}+\cdots +{\lambda }_{b-2}^{-1}{{\Pi}}_{b-1,b-1}+{\lambda }_{b-1}^{-1}}}.\hfill \end{align} \tag{ 19 }$$

For a process that starts at a transient initial state m with 0 < m < b, i.e. away from the reflecting state, the lower bound of the first sum in equation (19) has to be changed from i = 0 to i = m. Explicitly writing out the double sum in this way shows that terms of the outer sum (summation over i) correspond to single MFPTs ${\mathcal{T}}_{1}(i;0,i+1)$ from state i to i + 1 while keeping the reflective state unaltered in state 0. That is, irrespective of where the systems starts, it is always allowed to visit state '0' again. For illustration, consider the single MFPT ${\mathcal{T}}_{1}(2;0,3)$ for the transition 2 → 3, i.e. the case i = 2 above. The last contribution, ${\lambda }_{2}^{-1}$, is the time needed for the direct path from state 2 to 3. Depending on the transitions rates (e.g. when μ₂ > λ₂) the system might however also transition backwards to state 1 or even 0 before finally reaching the absorbing target in state 3 for the first time. These additional paths are accounted for by the corresponding inverse rates and the weighting factor Π_2,2 as defined in equation (18). Consequently, we may express the MFPT as a sum over all individual single MFPTs as

$$\begin{equation}{\mathcal{T}}_{1}(0;0,b)=\sum\limits _{i=0}^{b-1}{\mathcal{T}}_{1}(i;0,i+1)={\mathcal{T}}_{1}(0;0,1)+{\mathcal{T}}_{1}(1;0,2)+\cdots +{\mathcal{T}}_{1}(b-1;0,b),\end{equation} \tag{ 20 }$$

where the systems starts initially in state 0. In the same way as before, if one starts in a state 0 < m < b the summation starts from i = m.

A decomposition into single MFPTs now makes it possible to distinguish between larger and smaller contributions to the full MFPT. Single MFPTs from state i to state i + 1 (leaving the reflecting state 0 unchanged) are shown in figures 7(c) and (d) for the network description of the NM2-dimer system for antiparallel and parallel configuration, respectively. It becomes immediately evident that the single MFPTs are highly irregular and vary over several orders of magnitude depending on the initial state i₀. Based on this coarse-grained description we are now able to identify the largest contributions towards the average detachment time of the NM2 dimer (marked in green). This implies that if we were to initially prepare the system in state m > i^⋆, where i^⋆ denotes the state with the largest single MFPT ${\mathcal{T}}_{1}(i;0,i+1)$, we would expect a substantial decrease in the time needed until the NM2-dimer detaches as all contributions (and therefore the largest too) of single MFPTs for states < m drop out, i.e., the sum in equation (20) starts at i = m. Analogously, the next and all further drops in the MFPT can be identified using the same principle, namely finding the next largest single MFPT and choosing the initial state accordingly. In figures 7(c) and (d) we show that indeed our procedure of identifying the relevant minima and barriers (green bars) predicts the observed drops in the full solutions (green arrows and dotted lines), both for parallel (c) and antiparallel (d) configurations. Thus our main question has now been answered, namely which energy barriers lead to the large drops in the MFPTs.