Unfolding force definition and the unified model for the mean unfolding force dependence on the loading rate

Two sets of BD simulations were carried out. In the first set, the quartic free energy profile was ramped with constant loading rates ranging from 10 pN s⁻¹ to 10¹² pN s⁻¹ to the final force, which was 15 pN for low loading rates (equation (11)). In the second set of simulations, the Morse free energy profile was pulled through a spring with a spring constant of 100 pN nm⁻¹ with a constant velocity ranging from 1000 nm s⁻¹ to 10¹⁴ nm s⁻¹ to the final force, which was 150 pN for the low loading rates (equation (12)). The friction coefficient was chosen to be 10⁻⁵ pN s nm⁻¹ [35] in both sets of simulations, and the overdamped Langevin equation (equation (13)) was solved numerically using Euler–Maruyama method.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G\left(x,t \right)=100\left({{\left(\frac{x}{5} \right)}^{4}}-{{\left(\frac{x}{5} \right)}^{2}} \right)+4x-\overset{.}{\mathop{F}}\,tx\nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G\left(x,t \right)=41{{\left(1-{{e}^{-2\left(\frac{x}{0.1}-1 \right)}} \right)}^{2}}-41+\frac{\kappa}{2}{{\left(Vt-x+0.1 \right)}^{2}}.\nonumber \end{align} \tag{ 12 }$$

In equations (11) and (12), $t$ is the time and $V$ is the speed with which spring is pulled.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \zeta \overset{.}{\mathop{x}}\,\left(t \right)=-\frac{\partial G\left(x,t \right)}{\partial x}+\sqrt{2\zeta {{k}_{{\rm B}}}T}R\left(t \right).\nonumber \end{align} \tag{ 13 }$$

In equation (13), $R\left(t \right)$ is a delta-correlated Gaussian white noise with 0 mean. The simulation time step was chosen to be 10⁻⁹ s for the quartic free energy profile and 10⁻¹² s for the Morse free energy profile, securing less than 10% of error due to finite time step [36]. For high loading rates, where the simulation time becomes quite short, the time steps were decreased even more to ensure that at least half a million steps were in a trajectory (based on random walk as a limit of Brownian motion, Donker's theorem). For each loading rate, 3000 trajectories were generated. For each trajectory, the initial value for the reaction coordinate was chosen randomly from the Boltzmann distribution for the corresponding free energy profile's native state neighborhood. The first, reversible, and last unfolding forces have been recorded for each trajectory and the mean unfolding forces were calculated from that data.

To compare how well the models predict the outcome of BD simulations, was necessary to use the same parameters in the models that were used in the BD simulation. For the BD simulations with the quartic free energy profile, most of the parameters could be determined straightforwardly. The determination of the more 'complicated' parameters is described here. The energy profile parameter (µ) for the profile assumption (2016) model [23] was determined by fitting the free energy profile model proposed in [23] to the part of the quartic free energy profile. The unfolding rate was determined using equation (14), figure 1 [21, 37–39].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle k=\frac{{{k}_{{\rm B}}}T}{\zeta}\cdot \frac{1}{\int_{{{x}_{\min}}}^{b}{{{e}^{\frac{G\left(y \right)}{{{k}_{{\rm B}}}T}}}}\int_{a}^{y}{{{e}^{-\frac{G\left(z \right)}{{{k}_{{\rm B}}}T}}}{\rm d}z{\rm d}y}}.\nonumber \end{align} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

In equation (14), $G\left({} \right)$ is the free energy profile and a and b are reflecting and absorbing boundaries, respectively. The refolding rate was calculated using a similar equation. The equilibrium force was determined by assuming the force dependence of the free energy profile in equation (14) and numerically solving the equation of the unfolding and refolding rates in order to find the force at which they are equal.

In the case of the Morse free energy profile, the determination of the parameters is less trivial since this profile does not have any barrier in the absence of the external force. In order to find the distance between the native and transition states for this model, we choose as a barrier position 0.262 nm, which is the mean of the barrier positions when it first appears at 0.135 nm and when it disappears at 0.389 nm. Since the native state is at 0.1 nm, then ${{x}^{{\ddagger}}}$  ≈  0.162 nm. The unfolded state distance was chosen at 0.6 nm, which is where the second minimum becomes noticeable. Therefore, for this case the first unfolding force was taken as the product of the time at which system visited the unfolded state (0.6 nm) for the first time and the loading rate, the last unfolding force was taken as the product of the time at which the system was in the native state (0.1 nm) for the last time and the loading rate, and the reversible unfolding force was taken as the product of the loading rate and the time the system spent in the native state's well (time spent for x  <  0.262 nm). To determine the distance between the unfolded state and the transition state $x_{{\rm ref}}^{{\ddagger}}$ and the free energy difference between the transition state and the unfolded state $\Delta G_{{\rm ref}}^{{\ddagger}}$, for the Morse free energy profile, we followed the approach described in [40]. To determine the unfolding and refolding rates we used equation (14) and figure 1. The equilibrium force was determined with the same method as in the case of the quartic free energy profile. The determination of the parameter values for the Morse free energy profile is challenging since some parameters are virtual. Our choice of the Morse free energy profile with a spring for this study is due to the historical reasons [32]. We could have chosen a different path for determining the parameters of the Morse free energy profile, particularly by evaluating the first and third derivatives at the inflection point as described in [15]. However, parameters determined with this approach might put profile assumption models and related models in a favorable position. Ideally ${\kappa V}/{Z}$ should be used to calculate the loading rate [4, 34]. However, we have used $\kappa V$ since in our case $Z\approx 1.003 05$ which means that the soft spring approximation is applicable.

As a goodness of prediction criterion, we used R² as defined by equation (15).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}^{2}}\equiv 1-\frac{\sum\nolimits_{i=1}^{n}{{{\left({{O}_{i}}-{{E}_{i}} \right)}^{2}}}}{\sum\nolimits_{i=1}^{n}{{{\left({{O}_{i}}-\left\langle O \right\rangle \right)}^{2}}}}.\nonumber \end{align} \tag{ 15 }$$

In equation (15), ${{O}_{i}}$ is the observed value of the ith data point, i.e. the value we get from the data, ${{E}_{i}}$ is the expected value for the ith data point, i.e. the value predicted by the model, and $\left\langle O \right\rangle$ is the mean of the data (see figure 2(a)). Clearly, the closer the R² value is to 1, the better the model predicts the data. A model with more parameters has an advantage when models are being fitted. In such a case the adjusted R² might be used, which is defined by equation (16).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\bar{R}}^{2}}\equiv 1-\left(1-{{R}^{2}} \right)\frac{n-1}{n-p-1}.\nonumber \end{align} \tag{ 16 }$$

Zoom In Zoom Out Reset image size

In equation (16), ${{R}^{2}}$ is defined by equation (15), $n$ is the sample size (number of data points), and $p$ is the number of parameters in the model.

Another criterion used in this study for the goodness of the prediction is χ² defined by equation (17).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\chi}^{2}}\equiv \sum\limits_{i=1}^{n}{\frac{{{\left({{O}_{i}}-{{E}_{i}} \right)}^{2}}}{\sigma _{i}^{2}}}.\nonumber \end{align} \tag{ 17 }$$

In equation (17), ${{\sigma}_{i}}$ is the standard deviation for each data point and the summation goes through all the data points from first to the last (nth) (see figure 2(b)). It is important to note that instead of the standard deviation, the standard error or other types of errors can also be used to define χ². However, since all our data points are obtained from the averaging of the same number of unfolding forces, we use the standard deviation. Clearly, the closer χ² is to 0, the better the model predicts the data. Similar to adjusted R², there is also a reduced χ² that takes into account the degrees of freedom of the model, as defined by equation (18).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\bar{\chi}}^{2}}\equiv \frac{{{\chi}^{2}}}{n-p-1}.\nonumber \end{align} \tag{ 18 }$$

In equation (18), ${{\chi}^{2}}$ is given by equation (17). In our analyses we have not used the adjusted R² or reduced χ² since we did not fit models to the data, i.e. we did not change the parameters of the models to optimize R² or χ². By using the values in tables 1 and 2 and by determining the number of the data points for each range from figures 5 and 6, the adjusted R² and reduced χ² can easily be calculated for all models discussed in this work using equations (16) and (18).

Table 1. To quantify how well the models predict the simulated data we calculate R² and χ² for each model plotted in figure 5 at three regimes separately, and then all together (last column). Some of the models predict complex values for the mean unfolding force at high loading rates. The '—' in the table means that R² and χ² have complex values. The perfect prediction would result in R²  =  1 and χ²  =  0.

	Low loading rates ($\overset{.}{\mathop{F}}\,$  <  5000 pN s−1)		Intermediate loading rates (5000 pN s−1 $\leqslant$ $\overset{.}{\mathop{F}}\,$ $\leqslant$ 106 pN s−1)		High loading rates ($\overset{.}{\mathop{F}}\,$  >  106 pN s−1)		All loading rates
Model	R2	χ2	R2	χ2	R2	χ2	R2	χ2
Phenomenological (1997)	0.9695	0.1809	0.0275	12.6468	−0.4443	1882.94	−0.1878	1895.77
Profile assumption (2006)	0.8446	6.3365	—	—	—	—	—	—
Profile assumption (2008)	0.9889	0.2245	—	—	—	—	—	—
Profile assumption (2016)	0.842	6.2111	—	—	—	—	—	—
Rapid (2014)	0.9968	0.0187	0.9727	0.3191	0.9997	4.145	0.9997	4.4828
Unified first (2019)	0.9987	0.0093	0.997	0.0693	0.9998	0.4946	0.9998	0.5732
Unified reversible (2019)	0.6972	0.5088	0.9977	0.0553	0.9998	0.4961	0.9998	1.0602
Unified last (2019)	0.5613	2.2316	0.9952	0.1886	0.9998	0.4974	0.9998	2.9176

Table 2. To quantify how well the models predict the simulated data we calculate R² and χ² for each model plotted in figure 6 at three regimes separately and then all together (last column). The profile assumption (2006, 2008, 2010) model predicts complex values for the mean unfolding force at high loading rates. The '—' in the table means that R² and χ² have complex values. Perfect prediction would result in R²  =  1 and χ²  =  0.

	Low loading rates ($\overset{.}{\mathop{F}}\,$  <  107 pN s−1)		Intermediate loading rates (107 pN s−1 $\leqslant$ $\overset{.}{\mathop{F}}\,$ $\leqslant$ 1010 pN s−1)		High loading rates ($\overset{.}{\mathop{F}}\,$  >  1010 pN s−1)		All loading rates
Model	R2	χ2	R2	χ2	R2	χ2	R2	χ2
Phenomenological (1997, 2000, 2012)	−0.8937	8.5561	0.7867	8.6606	−0.4415	2534.49	−0.1848	2550.6
Profile assumption (2006, 2008, 2010)	0.6993	2.5944	—	—	—	—	—	—
Rapid (2014)	0.6861	2.726	0.9396	1.2329	0.9995	4.8118	0.9996	8.459
Unified first (2019)	0.6993	2.5944	0.8855	3.0696	0.9996	1.2131	0.9997	6.5727
Reversible (2012, 2016)	0.3416	1.2726	0.6724	4.6319	−0.4426	2554.35	−0.1858	2560.25
Unified reversible (2019)	−0.6159	2.7631	0.8735	4.1891	0.9996	1.2143	0.9997	8.1665
Unified last (2019)	−3.0559	1.5944	0.8723	4.275	0.9996	1.2144	0.9997	7.8419

The determination of the unfolding force for the case of a single transition (figure 3) is straightforward. However, it is unclear what force should be taken as an unfolding force if the system hops back and forth between various states before the final unfolding (figure 4). For such situations we define three unfolding forces. The first unfolding force is the product of the loading rate and the time when the system visits the unfolded state for the first time. The last unfolding force is the product of the loading rate and the time when the system was in the native state for the last time. The reversible unfolding force is the product of the loading rate and the time that the system spent in the native state's well (see yellow stripes in figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen from figure 3, the positions of the extrema of the free energy profile change with time if the external, time-dependent force is applied. However, the native, transition, and unfolded states are specified as particular values of the reaction coordinate and are independent of the externally applied force. In the above-stated example the native, transition, and unfolded states coincide in the absence of the external force with the first minimum, maximum, and the second minimum of the free energy profile. The state is specified by a particular value (or range of values) of the reaction coordinate, and it does not depend on the external force. What depends on the external force is the likelihood of finding the system in a particular state.

In order to see how the means of these three unfolding forces depend on the loading rate, we carry out BD simulations using quartic free energy profile under the external force that increases with a constant loading rate. We generate 3000 trajectories for each loading rate, for the loading rates that span 12 orders of magnitude (see figure 5). For further details, see the methods section.

Zoom In Zoom Out Reset image size

The surprising observation is that the mean last unfolding force shows non-monotonic dependence on the loading rate. At sufficiently low loading rates, the mean last unfolding force starts to increase with the decrease of the loading rate, a phenomenon that was not observed previously. This means that the time when the system was in the native state for the last time increases faster than the loading rate decreases, as a result the product of two increases. The behavior of the mean last unfolding force at even lower loading rates will be discussed at the end of this section (see figure 8). The mean reversible unfolding force at the low loading rate regime with the decrease of the loading rate tends to the constant value of F_eq. This means that the mean time spent in the native state's well increases to the same extent as the loading rate decreases. The decrease in the mean first unfolding force at the low loading rate regime is evident. As a reminder, the first unfolding force is the time when the first unfolding occurred multiplied by the loading rate. In the absence of the external force, there is a constant mean time (mean first passage time) when the system visits the unfolded state for the first time. In the presence of the slowly changing external force, the time when the first unfolding occurred tends to the mean first passage time. Hence, one of the two terms in the product tends to the constant while the other one decreases. This results in the decrease of the first unfolding force with the decreasing loading rate at the near-equilibrium regime.

As can be seen in figure 5, three differently defined mean unfolding forces converge with the increase of the loading rate. This is the result of moving away from the equilibrium and entering the regime where there are mostly single transitions in trajectories. In the case of a single transition, the difference between three unfolding forces is negligible as long as the time spent during the transition (transition path time) is negligible compared to the time spent before the transition (residence time of the native state). This is true for low loading rates for which also the hopping is common. This is not true for very high loading rates. For these rates only the mean reversible unfolding force is shown as the vast majority of trajectories at high loading rates have a single transition. For the single transition the determination of the unfolding force is straightforward: it is the time spent in the native state's well multiplied by the loading rate.

The next point that can be noticed from figure 5(b) is the kink around 5 · 10⁶ pN s⁻¹ of the loading rate value. After the critical value of the loading rate, the data points lie on a line with a higher slope. This is the indication of the change in the regime of unfolding. At such high loading rates, the unfolding mostly occurs after the disappearance of the barrier.

Hereby, we have specified three regimes: the low loading rate regime (close to equilibrium), the intermediate loading rate regime (non-equilibrium), and the high loading rate regime (non-equilibrium, ballistic unfolding). The smoothness of the kinks between these regimes (figure 5(b)) is due to the fact that the mean unfolding force for the given loading rate contains data from many trajectories. These trajectories can show multiple transitions between states, or can have a single transition to the unfolded state before or after the disappearance of the barrier.

To accurately predict the dependence of the mean unfolding force on the loading rate in all the three regimes, we developed a unified model. We used the force-dependent unfolding rate expression proposed in [15], and, following the method proposed in [22], we calculated the mean first unfolding force for the low and intermediate loading rates. Next, we extended that result for the case of high loading rates, when the unfolding is happening after the disappearance of the barrier [14, 18]. This results in equation (19), the derivations of which can be found in appendix A.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm First}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-{{\left| 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{\ddagger}}}{{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}}}{{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,} \right) \right|}^{\nu}} \right)\cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}-\overset{.}{\mathop{F}}\, \right) \nonumber \\ &+ \left(\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}} \right). \end{array}\nonumber \end{align} \tag{ 19 }$$

In equation (19), ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}}} \right)$, ${{k}_{{\rm 0}}}$ is not an independent parameter, and it can be calculated using the Kramers' approximation if the free energy profile is known. For example, for the linear-cubic profile ($\nu$  =  2/3), we will have ${{k}_{{\rm 0}}}=\frac{3\Delta {{G}^{{\ddagger}}}}{\pi \zeta {{x}^{{\ddagger ^{2}}}}}{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}}}$ in contrast to equation (5), where the harmonic-cusp profile is assumed. Generally using the linear-cubic free energy profile is more reasonable as discussed in [15, 34], and we will assume that profile in the future, i.e. $\nu$  =  2/3, if not specified. Here $\theta \left({} \right)$ is the Heaviside step function: $\newcommand{\e}{{\rm e}} \theta \left(x \right)=\left\{\begin{matrix} 1, & x>0 \\ 0, & x\leqslant 0 \end{matrix} \right.$. Hereby, the dependence of the mean first unfolding force on the loading rate is determined by 3 independent parameters: ${{x}^{\ddagger}}$, $\Delta {{G}^{\ddagger}}$, and $\zeta$.

Following the same procedure, except this time calculating the mean differently [32], we obtained the mean reversible unfolding force: equation (20) (see appendix A for details).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-\left| {{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}{{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu} \right) \nonumber \\ & \cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}-\overset{.}{\mathop{F}}\, \right)+\left(\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}} \right). \end{array}\nonumber \end{align} \tag{ 20 }$$

In equation (20), ${{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \right)$. Hereby, the dependence of the mean reversible unfolding force on the loading rate is determined by 4 independent parameters: ${{x}^{{\ddagger}}}$, $\Delta {{G}^{{\ddagger}}}$, ${{F}_{{\rm eq}}}$, and $\zeta$.

Notice that the critical loading rate depends on the definition of the unfolding force. Generally, the mean reversible unfolding force riches the given force at the lower loading rate than the mean first unfolding force since by definition for the same trajectory the first unfolding force is not higher than the reversible unfolding force. With the increase of the loading rate, the number of trajectories with hopping decreases, and, due to this, the first and reversible mean unfolding forces are getting closer, but eventually the mean reversible unfolding force riches the critical force at a slightly smaller loading rate than the mean first unfolding force.

The mean last unfolding force can be determined based on the symmetry arguments. In figure 4 we observe that $\left\langle {{F}_{{\rm Last}}} \right\rangle =\left\langle {{F}_{{\rm Rev}}} \right\rangle +\left\langle {{F}_{x}} \right\rangle$ where $\left\langle {{F}_{x}} \right\rangle$ is the mean force due to the time spent in the unfolded state's well (above 0.51 nm), before the last unfolding. $\left\langle {{F}_{x}} \right\rangle$ equals the time spent in the unfolded state's well before the last unfolding multiplied by the loading rate. Interestingly, at close to equilibrium regime, it is possible to calculate $\left\langle {{F}_{x}} \right\rangle$ from the relaxation protocol, namely, $\left\langle {{F}_{x}} \right\rangle \approx \left\langle {{F}_{{\rm FirstRef}}} \right\rangle -\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ where $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ is the man first refolding force, and $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ is the mean reversible refolding force. Hence, the mean last unfolding force can be calculated through equation (21):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{F}_{{\rm Last}}} \right\rangle \approx \left\langle {{F}_{{\rm Rev}}} \right\rangle +\left\langle {{F}_{{\rm FirstRef}}} \right\rangle -\left\langle {{F}_{{\rm RevRef}}} \right\rangle .\nonumber \end{align} \tag{ 21 }$$

At the non-equilibrium regime, the symmetry will break, but the equation (21) will still hold since in the case of a single transition the last two terms will cancel each other. $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ and $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ can be calculated following the same path as for $\left\langle {{F}_{{\rm First}}} \right\rangle$ and $\left\langle {{F}_{{\rm Rev}}} \right\rangle$ only using the force dependent refolding rate as in [41] and different integration limits (see appendix A for details). This results in equations (22) and (23):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left(\left| {{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}{{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu}-1 \right) \nonumber \\ & \cdot\theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}} \right). \end{array}\nonumber \end{align} \tag{ 22 }$$

In equation (22), ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right)$, ${{k}_{{\rm 0ref}}}=\frac{3\Delta G_{{\rm ref}}^{{\ddagger}}}{\pi \zeta x{{_{{\rm ref}}^{{\ddagger ^{2}}}}}}{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}}}$, $x_{{\rm ref}}^{{\ddagger}}$ is the distance between the unfolded state and the transition state, ${{k}_{{\rm 0ref}}}$ is the refolding rate in the absence of the external force, $\Delta G_{{\rm ref}}^{{\ddagger}}$ is the free energy difference between the transition state and the unfolded state, and ${{F}_{{\rm fin}}}$ is the final force till which the system is pulled and from which the (hypothetical) refolding protocol starts. The final force is set by the experimental design and is not an intrinsic parameter describing the system.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm RevRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left(\left| {{\left(1+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}{{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu}-1 \right)\cdot \nonumber \\ & \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}} \right). \end{array}\nonumber \end{align} \tag{ 23 }$$

In equation (23), ${{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(1+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right)$. In this case, ${{F}_{{\rm eq}}}$ is not an independent parameter since in equations (22) and (23) there are parameters involved which describing the unfolded state's well. The equilibrium unfolding force can be calculated using the probability of having the system folded at a given force [33]: ${{F}_{{\rm eq}}}=\int_{0}^{\infty}{\frac{1}{1+{{e}^{\frac{F\left({{x}^{{\ddagger}}}+x_{{\rm ref}}^{{\ddagger}} \right)-\left(\Delta {{G}^{{\ddagger}}}-\Delta G_{{\rm ref}}^{{\ddagger}} \right)}{{{k}_{{\rm B}}}T}}}}}=\frac{{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}+x_{{\rm ref}}^{{\ddagger}}}\ln \left(1+{{e}^{\frac{\Delta {{G}^{{\ddagger}}}-\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}}} \right)$. This formula of ${{F}_{{\rm eq}}}$ is model independent in the sense that we do not assume force dependent unfolding and refolding rates. By inserting equations (20), (22) and (23) into (21), the dependence of the mean last unfolding force on the loading rate is determined by 5 independent parameters: ${{x}^{{\ddagger}}}$, $\Delta {{G}^{{\ddagger}}}$, $x_{{\rm ref}}^{{\ddagger}}$, $\Delta G_{{\rm ref}}^{{\ddagger}}$, and $\zeta$. Independently of applying or not applying relaxation protocol, equation (21) (which includes equation (20), (22), and (23)) is applicable for predicting the dependence of the mean last unfolding force on the loading rate.

Exponential integral and Heaviside step function can be approximated using other well-known functions, for example, ${{e}^{x}}{{E}_{1}}\left(x \right)\approx \ln \left(1+\frac{{{e}^{-\gamma}}}{x} \right)$ and $\theta \left(x \right)\approx \frac{1}{1+{{e}^{-qx}}}$ where q is a large positive number [11]. However, for maximum precision, the fitting function files provided in the supplementary material (available online at stacks.iop.org/JSTAT/2020/033201/mmedia) exclude these approximations. These files can be directly used for data analysis.

To show how well all above-described models predict the dependence of the mean unfolding force on the loading rate, the models' predictions and the data from the BD simulations are presented on the same graph (figure 5). The parameters in these models are the same as in the BD simulations. To further quantify how well the models predicted the BD simulated data, R² and χ² were calculated (see table 1). For further details, see the methods section.

We observe that the phenomenological model [5, 9] fails at high and intermediate loading rates. The main disadvantage of the profile assumption models [15, 22, 23] is that they fail at high loading rates. The rapid model [18], which combines the result of the phenomenological model [9] with the high loading rate result [14, 18], arguably provides the best prediction among the previously proposed models. The reason for this could be that the intermediate loading rate regime range is small compared to the high loading rate regime, and more data is generated in the high loading rate regime. The loading rate range for the BD simulations is chosen based on the previous experimental and molecular dynamics simulations studies [13]. The unified first model clearly gives the overall best prediction. According to unified last model, the mean last unfolding force tends to the constant at very low loading rates as it will be shown later in figure 8. However, conducting simulations at very low loading rates is challenging due to long simulation times and large memory required to deal with simulated data. Hence, in the case of very low loading rates, the prediction of the unified last model needs further verification.

We further extend the unified model for the case when the system is being pulled through the spring. The dependence of the mean unfolding forces on the loading rate is derived using the same methods as when the force is applied directly. The only difference here is that different force dependent unfolding and refolding rates are used (see appendix B for derivations) [34]. The calculations result in equations (24)–(27) for the mean first unfolding force, mean reversible unfolding force, mean first refolding force, and mean reversible refolding force respectively. The last three are used for calculating the mean last unfolding force in equation (21).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm First}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-{{\left| {{Z}^{\frac{2}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}{{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{Z}^{\frac{2}{\nu}}} \right)}}}}{{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{Z}^{\frac{2}{\nu}}} \right)}} \right) \right|}^{\nu}} \right) \nonumber \\ &\cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}-\overset{.}{\mathop{F}}\, \right) + \left(\frac{\Delta {{G}^{{\ddagger}}}Z}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}} \right). \end{array}\nonumber \end{align} \tag{ 24 }$$

In equation (24), ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}Z}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}{{Z}^{\frac{2}{\nu}}}}{{{k}_{{\rm B}}}T}}} \right)$, and Z is given by equation (9). For the linear-cubic free energy profile $Z=1+\frac{\kappa {{x}^{{\ddagger ^{2}}}}}{6\Delta {{G}^{{\ddagger}}}}$.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-\left| {{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}{{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu} \right) \nonumber \\ & \cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}-\overset{.}{\mathop{F}}\, \right)+\left(\frac{\Delta {{G}^{{\ddagger}}}Z}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}} \right) \end{array}\nonumber \end{align} \tag{ 25 }$$

In equation (25), ${{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}Z}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \right)$.

$\left\langle {{F}_{{\rm Last}}} \right\rangle$ could be calculated using equations (21), (25), (26) and (27) for $\left\langle {{F}_{{\rm Rev}}} \right\rangle$, $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ and $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ respectively.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(\left| {{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}{{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu}-Z_{{\rm ref}}^{2} \right)\cdot \nonumber \\ & \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}} \right). \end{array}\nonumber \end{align} \tag{ 26 }$$

In equation (26), ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right)$, and ${{Z}_{{\rm ref}}}=1+\frac{\kappa}{{{\kappa}_{{{{\rm G}}_{{\rm 0ref}}}}}\left({{x}_{\min}} \right)}$ where $\newcommand{\e}{{\rm e}} {{\kappa}_{{{{\rm G}}_{{\rm 0ref}}}}}\left({{x}_{\min}} \right)\equiv {{\left. \frac{{{\partial}^{2}}{{G}_{{\rm 0ref}}}\left(x \right)}{\partial {{x}^{2}}} \right|}_{x={{x}_{\min}}}}$ is the curvature of the free energy profile in the unfolded state, and ${{G}_{{\rm 0ref}}}\left(x \right)=\Delta G_{{\rm ref}}^{^{{\ddagger}}}\frac{\nu}{1-\nu}{{\left(\frac{x}{x_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{1-\nu}}}-\frac{\Delta G_{{\rm ref}}^{^{{\ddagger}}}}{\nu}\frac{x}{x_{{\rm ref}}^{{\ddagger}}}$ is the model for the free energy profile. For the linear-cubic free energy profile ${{Z}_{{\rm ref}}}=1-\frac{\kappa x_{{\rm ref}}^{{\rm 2}}}{6\Delta G_{{\rm ref}}^{{\ddagger}}}$.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm RevRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(\left| {{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}{{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \nonumber \\ &\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \right|^{\nu}-Z_{{\rm ref}}^{2} \right)\cdot \nonumber \\ & \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}} \right). \end{array}\nonumber \end{align} \tag{ 27 }$$

In equation (27), ${{\overset{.}{\mathop{F}}\,}_{{\rm cRevRef}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm eq}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right)$. As previously ${{F}_{{\rm eq}}}$ (in the case of the mean last unfolding force) is not an independent parameter and can be expressed through other parameters of the system [33]

${{F}_{{\rm eq}}}=\int_{0}^{\infty}{\frac{1}{1+{{e}^{\frac{\frac{{{F}^{2}}}{2\kappa}-\left(\Delta {{G}^{{\ddagger}}}-\Delta G_{{\rm ref}}^{{\ddagger}} \right)}{{{k}_{{\rm B}}}T}}}}}=-\sqrt{\frac{\pi \kappa {{k}_{{\rm B}}}T}{2}}{\rm L}{{{\rm i}}_{{1}/{2}}}\left(-{{e}^{\frac{\Delta {{G}^{{\ddagger}}}-\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}}} \right)$ where ${\rm L}{{{\rm i}}_{{1}/{2}}}\left({} \right)$ is the polylogarithm [11].

In order to see how well different models predict the dependence of the mean unfolding force on the loading rate when a pulling force is applied through a spring, we carried out BD simulations using Morse free energy profile which was pulled with constant speed through a spring. The spring constant was 100 pN nm⁻¹ and the Morse free energy profile parameters were chosen as in [32] (see methods section for details).

To show how well above discussed models predict the dependence of the mean unfolding force on the loading rate, we plot the models' predictions and the data from the BD simulations on the same graph (figure 6). The parameters for plotting the models were the same as those used in the BD simulations. To quantify how well the models predict the BD simulation data, we have calculated R² and χ² (see methods section for details), which are presented in table 2. Unified models generally give the best predictions. It appears that the soft spring approximation works well in this case because Z is close to 1, and hence we could have used models that do not take into account the correction for the spring. However, since the reversible model [32, 33] was originally developed assuming pulling through the spring, we used models that take into account the effect of the spring in order to have all models on equal ground.

Zoom In Zoom Out Reset image size

The distributions of differently defined unfolding forces in the low loading rate regime are vastly different (figure 7). The reversible unfolding force distribution is symmetric and narrow. The last unfolding force is more widely distributed and has a right tail. Generally at intermediate loading rates, in cases of single transition, the unfolding force distribution has a left tail [14–18]. When there is hopping, the first unfolding force distribution also shows a left tail, but as the loading rate decreases the overall distribution moves closer to 0 and the left tail disappears, as can be seen in figure 7(a). The evolution of the differently defined unfolding force distributions over the loading rate is an interesting question since the distributions are very different in the near-equilibrium regime but eventually coincide as the system moves further from equilibrium. A subject of future research could be the derivation of a differential equation describing these 'dynamics'.

Zoom In Zoom Out Reset image size

As described above, in the low loading rate regime the mean last unfolding force starts to increase with the decrease of the loading rate. Performing BD simulations at lower loading rates is challenging because it requires long simulation times and large memory. However, we could see the prediction of the unified last model at lower loading rates (figure 8).

Zoom In Zoom Out Reset image size

We first observe an increase and then a plateau of the mean last unfolding force at low loading rates. This means that the time for the last unfolding increases with decreasing loading rate, such that the product of that time and the loading rate still increases. When the loading rate goes to 0, $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ (equation (22)) tends toward ${{F}_{{\rm fin}}}$, i.e. to the final force until which the system is pulled or the force from which the relaxation could start. At the same limit, $\left\langle {{F}_{{\rm Rev}}} \right\rangle$ (equation (20)) and $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ (equation (23)) both tend toward the same constant ${{F}_{{\rm eq}}}$, hence $\left\langle {{F}_{{\rm Last}}} \right\rangle$ (equation (21)) will tend toward ${{F}_{{\rm fin}}}$ as $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$. The plateau is an indication that at a sufficiently low loading rates hopping will continue until the end of the process, i.e. until the final force to which the system is driven. This occurs because there is always a finite probability that the system will refold back and unfold again, even if the energy of the native state has become much higher than that of unfolded state during the pulling process. Over sufficiently long time periods, such as those needed to accomplish the entire pulling process at very low loading rates, such improbable events are likely to occur.

The derivation of equation (20) will be shown in two steps. First, the dependence of the mean unfolding force on the loading rate for the low and intermediate loading rates is calculated. Second, that result is united with the result for high loading rates.

For the first part, we follow the approach of [32] with a difference that we use the profile assumption rate expression [15] (equation (A.1)) which is more general than the phenomenological rate [6–8] used in [32] and reduces to it when $\nu$  =  1.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle k\left(F \right)={{k}_{{\rm 0}}}{{\left(1-\frac{\nu F{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}-1}}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu F{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}\nonumber \end{align} \tag{ A.1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int{k\left(F \right){\rm d}F=\frac{{{k}_{0}}{{k}_{{\rm B}}}T}{{{x}^{\ddagger}}}}{{e}^{\frac{\Delta {{G}^{\ddagger}}}{{{k}_{B}}T}\left(1-{{\left(1-\frac{\nu F{{x}^{\ddagger}}}{\Delta {{G}^{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}.\nonumber \end{align} \tag{ A.2 }$$

Using the approach of [32], by substituting equation (A.2) into equation (A.3) we get equation (A.4).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{1}^{\,p}{\frac{{\rm d}{{p}^{\prime}}}{{{p}^{\prime}}}}=-\frac{1}{\overset{.}{\mathop{F}}\,}\int_{{{F}_{{\rm eq}}}}^{F}{k\left({{F}^{\prime}} \right)}{\rm d}{{F}^{\prime}}\nonumber \end{align} \tag{ A.3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln p=-\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}\left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu F{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right).\nonumber \end{align} \tag{ A.4 }$$

By solving equation (A.4) for F we find equation (A.5).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F=\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-{{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,\ln p}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}} \right).\nonumber \end{align} \tag{ A.5 }$$

The advantage of using the rate in equation (A.1) [15] versus the rate used in [23], where the first power is $2-1/\mu$ instead of $1/\nu -1$, is that the indefinite integral in equation (A.2) will contain an En-function [10, 11] which complicates the later step of solving equation (A.4) for F.

At this step, we follow the approximate averaging procedure as in [22].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \alpha \equiv \ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,\ln p}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right) \nonumber \\ & \left\langle F \right\rangle \approx F\left(\left\langle \alpha \right\rangle \right) \nonumber \\ & \int_{0}^{1}{\ln \left(A-\frac{\ln p}{B} \right)}{\rm d}p={{e}^{AB}}{{E}_{1}}\left(AB \right)+\ln \left(A \right)\approx \ln \left(A+\frac{{{e}^{-\gamma}}}{B} \right). \end{array}\nonumber \end{align} \tag{ A.6 }$$

From equations (A.5) and (A.6), we obtain equation (A.7) for the mean reversible unfolding force.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\left({{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right)^{\nu} \right).\nonumber \end{align} \tag{ A.7 }$$

Using the approximation from equation (A.6) in equation (A.7), we get equation (A.8).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-{{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,{{e}^{-\gamma}}}{{{k}_{{\rm 0un}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}} \right). \nonumber \end{align} \tag{ A.8 }$$

Note that equations (A.7) and (A.8) are valid until the disappearance of the barrier. Mathematically, this is demonstrated by the negative values taken to the non-integer power. Next, we determine the minimum loading rate at which this occurs.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,{{e}^{-\gamma}}}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right)=0 \nonumber \\ & {{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,{{e}^{-\gamma}}}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}={{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}}} \nonumber \\ & \overset{.}{\mathop{F}} = \frac{{{k}_{{\rm 0}}}{{k}_{B}}T{{e}^{\gamma}}}{{{x}^{{\ddagger}}}}\left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}}}-{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right) \nonumber \\ & \overset{.}{\mathop{F}} = \frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \right). \end{array}\nonumber \end{align}$$

At the critical loading rate, the mean reversible unfolding force is approximately $\Delta {{G}^{{\ddagger}}}/\nu {{x}^{{\ddagger}}}$ according to equation (A.7) and equation (A.8). After the barrier disappearance, the mean unfolding force increases as the square root of the loading rate $\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta}$ [14, 18]. In order to combine the two functions that are valid for the different loading rate regimes, we take the absolute value of the expression that is raised to the non-integer power in equation (A.7) and use the Heaviside step function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-\left| 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\left({{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right|^{\nu} \right) \nonumber \\ & \cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}-\overset{.}{\mathop{F}}\, \right)+\left(\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}} \right) \end{array}\nonumber \end{align} \tag{ A.9 }$$

where ${{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}} \bigg(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \bigg)$ and ${{k}_{{\rm 0}}}=\frac{3\Delta {{G}^{{\ddagger}}}}{\pi \zeta {{x}^{{\ddagger ^{2}}}}}{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}}}$.

Equation (A.9) can be further simplified to equation (20), or using the approximation from equation (A.6) to equation (A.10).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}\left(1-{{\left| 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{e}^{-\gamma}}{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,}{{{k}_{{\rm 0un}}}{{k}_{{\rm B}}}T} \right) \right|}^{\nu}} \right)\cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}-\overset{.}{\mathop{F}}\, \right) \nonumber \\ &+ \left(\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}} \right). \end{array}\nonumber \end{align} \tag{ A.10 }$$

The dependence of the mean first unfolding force on the loading rate equation (19) can be obtained by starting the integration in equation (A.3) from 0 and following the same steps, or by directly setting ${{F}_{{\rm eq}}}=0$ in equation (A.9) and further simplifying.

The mean first refolding force equation (22) and the mean reversible refolding force equation (23) that are used for calculating the mean last unfolding force in equation (21) can be obtained using the same methods as described above, though by using equation (A.11) for the refolding rate [41] and by using limits of integration for $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ from ${{F}_{{\rm fin}}}$ to F in equation (A.13) and for $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ from ${{F}_{{\rm eq}}}$ to $F$ in a similar equation.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{k}_{{\rm ref}}}\left(F \right)={{k}_{{\rm 0ref}}}{{\left(1+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}-1}}{{e}^{\frac{\Delta G_{{\rm ref}}^{^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}\nonumber \end{align} \tag{ A.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int{{{k}_{{\rm ref}}}\left(F \right)dF=-\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}.\nonumber \end{align} \tag{ A.12 }$$

By substituting equation (A.12) into equation (A.13) we get equation (A.14).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int\limits_{1}^{\,p}{\frac{{\rm d}{{p}^{\prime}}}{{{p}^{\prime}}}}=\frac{1}{\overset{.}{\mathop{F}}\,}\int_{{{F}_{{\rm fin}}}}^{F}{{{k}_{{\rm ref}}}\left({{F}^{\prime}} \right)}{\rm d}{{F}^{\prime}}\nonumber \end{align} \tag{ A.13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln p=-\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}\left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-{{e}^{\frac{\Delta G_{\rm ref}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right).\nonumber \end{align} \tag{ A.14 }$$

By solving equation (A.14) for F we find equation (A.15).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F=\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left({{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,\ln p}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}}-1 \right).\nonumber \end{align} \tag{ A.15 }$$

At this step, we follow the approximate averaging procedure as in [22].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha \equiv \ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,\ln p}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right)\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle F \right\rangle \approx F\left(\left\langle \alpha \right\rangle \right).\nonumber \end{align}$$

From equations (A.15) and (A.6), we obtain equation (A.16) for the mean first refolding force.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left(\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\left({{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right)^{\nu}-1 \right).\nonumber \end{align} \tag{ A.16 }$$

Using the approximation from equation (A.6) in equation (A.16), we get equation (A.17).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left({{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{e}^{-\gamma}}x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}}-1 \right).\nonumber \end{align} \tag{ A.17 }$$

Next, we find the critical loading rate in the same way as previously.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{e}^{-\gamma}}x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right)=0\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overset{.}{\mathop{F}} = \frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right).\nonumber \end{align}$$

By combining the two functions that are valid for the different loading rate regimes as previously, using the Heaviside step function we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}\left(\left| 1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\left({{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right|^{\nu}-1 \right)\cdot \nonumber \\ & \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}} \right) \end{array}\nonumber \end{align} \tag{ A.18 }$$

where ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirstRef}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(1+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right)$ and ${{k}_{{\rm 0ref}}}=\frac{3\Delta G_{{\rm ref}}^{{\ddagger}}}{\pi \zeta x{{_{{\rm ref}}^{{\ddagger ^{2}}}}}}{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}}}$.

Equation (A.18) can be further simplified to equation (22), or by using the approximation from equation (A.6) in a similar equation to equation (A.10).

The dependence of the mean reversible refolding force on the loading rate equation (23) can be obtained by starting the integration in equation (A.13) from ${{F}_{{\rm eq}}}$ and following the same steps, or by directly setting ${{F}_{{\rm fin}}}$ to ${{F}_{{\rm eq}}}$ in equation (A.18) and further simplifying.

When the force is being applied through the spring, we need to use an appropriate expression for the force-dependent unfolding rate equation (B.1) [34], which then can be plugged into equation (B.2).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle k\left(F \right)={{k}_{{\rm 0}}}{{\left({{Z}^{2}}-\frac{\nu F{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}-1}}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu F{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}\nonumber \end{align} \tag{ B.1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int{k\left(F \right){\rm d}F=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}Z}}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu F{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}.\nonumber \end{align} \tag{ B.2 }$$

Using the approach of [32], by substituting equation (B.2) into equation (B.3) we get equation (B.4).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{1}^{\,p}{\frac{{\rm d}{{p}^{\prime}}}{{{p}^{\prime}}}}=-\frac{1}{\overset{.}{\mathop{F}}\,}\int_{{{F}_{{\rm eq}}}}^{F}{k\left({{F}^{\prime}} \right)}{\rm d}{{F}^{\prime}}\nonumber \end{align} \tag{ B.3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln p=-\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}Z\overset{.}{\mathop{F}}\,}\left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu F{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right).\nonumber \end{align} \tag{ B.4 }$$

By solving equation (B.4) for F, we find equation (B.5).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F=\frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-{{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{{{x}^{{\ddagger}}}\overset{\, .}{\mathop{F}}Z\,\ln p}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}} \right).\nonumber \end{align} \tag{ B.5 }$$

At this step, we follow the approximate averaging procedure as in [22].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha \equiv \ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z\ln p}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right)\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle F \right\rangle \approx F\left(\left\langle \alpha \right\rangle \right).\nonumber \end{align}$$

From equations (B.5) and (A.6) we obtain the mean reversible unfolding force:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\left({{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}}\right.\right. \right.\nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right)^{\nu} \right).\nonumber \end{align} \tag{ B.6 }$$

Using the approximation from equation (A.6) in equation (B.6), we get equation (B.7).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-{{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z{{e}^{-\gamma}}}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}} \right).\nonumber \end{align} \tag{ B.7 }$$

Next, we find the critical loading rate in the same way as in appendix A.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\ln \left({{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z{{e}^{-\gamma}}}{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T} \right)=0\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overset{.}{\mathop{F}} = \frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}Z}\left(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \right).\nonumber \end{align}$$

By combining the two functions that are valid for the different loading rate regimes as was done in appendix A, using the Heaviside step function we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm Rev}}} \right\rangle \approx \frac{\Delta {{G}^{{\ddagger}}}}{\nu {{x}^{{\ddagger}}}Z}\left({{Z}^{2}}-\left| 1-\frac{{{k}_{{\rm B}}}T}{\Delta {{G}^{{\ddagger}}}}\left({{e}^{\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T}{{{x}^{{\ddagger}}}\overset{.}{\mathop{F}}\,Z}{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right|^{\nu} \right)\cdot \nonumber \\ & \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}-\overset{.}{\mathop{F}}\, \right)+\left(\frac{\Delta {{G}^{{\ddagger}}}Z}{\nu {{x}^{{\ddagger}}}}-\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}{{x}^{{\ddagger}}}\zeta}+\sqrt{2\overset{.}{\mathop{F}}\,{{x}^{{\ddagger}}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}} \right) \end{array}\nonumber \end{align} \tag{ B.8 }$$

where ${{\overset{.}{\mathop{F}}\,}_{{\rm cRev}}}=\frac{{{k}_{{\rm 0}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}+\gamma}}}{{{x}^{{\ddagger}}}Z} \bigg(1-{{e}^{-\frac{\Delta {{G}^{{\ddagger}}}}{{{k}_{{\rm B}}}T}{{\left({{Z}^{2}}-\frac{\nu {{F}_{{\rm eq}}}{{x}^{{\ddagger}}}Z}{\Delta {{G}^{{\ddagger}}}} \right)}^{\frac{1}{\nu}}}}} \bigg)$.

Equation (B.8) can be further simplified to equation (25), or by using the approximation from equation (A.6) to a similar equation to equation (A.10).

The dependence of the mean first unfolding force on the loading rate equation (24) can be obtained by starting the integration in equation (B.3) from 0 and following the same steps, or by directly setting ${{F}_{{\rm eq}}}=0$ in equation (B.8) and further simplifying.

The mean first refolding force equation (26) and the mean reversible refolding force equation (27) that are used for calculating the mean last unfolding force when the force is applied through the spring in equation (21) can be obtained using the same methods as described above, using equation (B.9) for the refolding rate [41] and limits of integration for $\left\langle {{F}_{{\rm FirstRef}}} \right\rangle$ from ${{F}_{{\rm fin}}}$ to $F$ in equation (B.11) and for $\left\langle {{F}_{{\rm RevRef}}} \right\rangle$ from ${{F}_{{\rm eq}}}$ to $F$ in a similar equation.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{k}_{{\rm ref}}}\left(F \right)={{k}_{{\rm 0ref}}}{{\left(Z_{{\rm ref}}^{2}+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}-1}}{{e}^{\frac{\Delta G_{\rm ref}^{^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}\nonumber \end{align} \tag{ B.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int{{{k}_{{\rm ref}}}\left(F \right){\rm d}F=-\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{^{{\ddagger}}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}.\nonumber \end{align} \tag{ B.10 }$$

By substituting equation (B.10) into equation (B.11) we get equation (B.12).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{1}^{\,p}{\frac{{\rm d}{{p}^{\prime}}}{{{p}^{\prime}}}}=\frac{1}{\overset{.}{\mathop{F}}\,}\int_{{{F}_{{\rm fin}}}}^{F}{{{k}_{{\rm ref}}}\left({{F}^{\prime}} \right)}{\rm d}{{F}^{\prime}}\nonumber \end{align} \tag{ B.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln p=-\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}\left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu Fx_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right).\nonumber \end{align} \tag{ B.12 }$$

By solving equation (B.12) for F, we find equation (B.13).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F=\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left({{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}\ln p}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}}-Z_{{\rm ref}}^{2} \right). \nonumber \end{align} \tag{ B.13 }$$

At this step, we follow the approximate averaging procedure as in [22].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha \equiv \ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}-\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}\ln p}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right)\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle F \right\rangle \approx F\left(\left\langle \alpha \right\rangle \right).\nonumber \end{align}$$

From equations (B.13) and (A.6), we obtain for the mean first refolding force equation (B.14).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\left({{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right. \right.\nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right)^{\nu}-Z_{{\rm ref}}^{2} \right).\nonumber \end{align} \tag{ B.14 }$$

Using the approximation from equations (A.6) in (B.14), we get equation (B.15).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left({{\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}{{e}^{-\gamma}}}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right) \right)}^{\nu}}-Z_{{\rm ref}}^{2} \right).\nonumber \end{align} \tag{ B.15 }$$

Next, we find the critical loading rate in the same way as previously.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\ln \left({{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}+\frac{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}{{e}^{-\gamma}}}{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T} \right)=0\nonumber \end{align}$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overset{.}{\mathop{F}} = \frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \right).\nonumber \end{align}$$

By combining the two functions that are valid for the different loading rate regimes as previously, using the Heaviside step function we obtain equation (B.16):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \left\langle {{F}_{{\rm FirstRef}}} \right\rangle \approx \frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{\nu x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}\left(\left(1-\frac{{{k}_{{\rm B}}}T}{\Delta G_{{\rm ref}}^{{\ddagger}}}\left({{e}^{\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}}}} \right.\right.\right. \nonumber \\ &\left.\left.\left. {{E}_{1}}\left(\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T}{x_{{\rm ref}}^{{\ddagger}}\overset{.}{\mathop{F}}\,{{Z}_{{\rm ref}}}}{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right)}} \right)+\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}\left(1-{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}} \right) \right) \right)^{\nu}-Z_{{\rm ref}}^{2} \right) \nonumber \\ &\cdot \theta \left({{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}-\overset{.}{\mathop{F}}\, \right)+\left(-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\nu x_{{\rm ref}}^{{\ddagger}}}+\sqrt{2{{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}x_{{\rm ref}}^{{\ddagger}}\zeta}-\sqrt{2\overset{.}{\mathop{F}}\,x_{{\rm ref}}^{{\ddagger}}\zeta} \right)\cdot \theta \left(\overset{.}{\mathop{F}} - {{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}} \right) \end{array}\nonumber \end{align} \tag{ B.16 }$$

where ${{\overset{.}{\mathop{F}}\,}_{{\rm cFirst}}}=\frac{{{k}_{{\rm 0ref}}}{{k}_{{\rm B}}}T{{e}^{\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}+\gamma}}}{x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}} \bigg (1-{{e}^{-\frac{\Delta G_{{\rm ref}}^{{\ddagger}}}{{{k}_{{\rm B}}}T}{{\left(Z_{{\rm ref}}^{2}+\frac{\nu {{F}_{{\rm fin}}}x_{{\rm ref}}^{{\ddagger}}{{Z}_{{\rm ref}}}}{\Delta G_{{\rm ref}}^{{\ddagger}}} \right)}^{\frac{1}{\nu}}}}} \bigg)$.

Equation (B.16) can be further simplified to equation (26) or using the approximation from equation (A.6) to a similar equation to equation (A.10).

The dependence of the mean reversible refolding force on the loading rate when the force is applied through a spring, given by equation (27), can be obtained by starting the integration in equation (B.11) from ${{F}_{{\rm eq}}}$ and following the same steps, or by directly setting ${{F}_{{\rm fin}}}$ to ${{F}_{{\rm eq}}}$ in equation (B.16) and further simplifying.