The time-changed stochastic approach and fractionally integrated processes to model the actin-myosin interaction and dwell times

1.   Introduction

1.1.   The phenomenon

The interaction between myosin and actin proteins inside the sarcomere is the basis of the muscle contraction process. The resultant mechanical work is essentially produced by the chemical energy released by ATP hydrolysis during successive Steps of the sliding of the myosin filaments. Usually, in the modeling context, we say that the myosin moves along the actin filament performing Steps mostly in a forward direction but sometimes also in a backward direction. The mechanism at the origin of this energy production is still not completely understood (see, for instance, ^[1,2,3] and references therein). More specifically, a skeletal muscle is composed by basic units, each called a sarcomere. Each sarcomere is essentially composed by two proteins: the actin and the myosin. The (deterministic) lever-arm theory (^[4]) states that the reciprocal sliding of these proteins leads to the contraction and relaxation of muscles. The so-called subfragment-1 (S1) is the structure of the head, neck, and light chains of myosin. It is reasonable to interpret this as an independent generator of force and movement (^[5,6]). Indeed, it is well-know that the molecular motor S1 acts as an enzyme that catalyzes the hydrolysis of ATP. This phenomenon occurs through a chemical and mechanical cycle: S1 attaches to the actin in the binding site, producing a cross-bridge, which contributes to hydrolyze ATP and causes a consequent conversion of chemical energy into movements and heat. Indeed, when the myosin head is bound to the actin, inorganic phosphate is released and the consequent power stroke takes place: this is the moment in which the myosin flexes pulling the actin along with it. The occurrence of the power stroke requires hydrolysis of ATP, which breaks a high-energy phosphate bond to release energy. The neck swing and the thermal bath in which the myosin is embedded determine the position of the myosin head and its potential energy, during the attachment of the myosin head to the actin filament and until to the final release of the phosphoric anion.

1.2.   Biological motivations and mathematical justifications

From a mathematical and modeling point of view, we aim to describe the dynamics of the myosin head in a binding site of the actin, during a cross-bridge, by means of the space-time evolution of a point-wise particle in a potential well; the model is a stochastic process $X(t)$ representing the position of the particle at time $t$ in a potential well. The potential of $X(t)$, denoted by $U(x)$, is a space-periodic function, and $U^\prime(x)$ stands for its potential energy. In the past, several potential profiles were investigated (^[7]). Here, we consider the parabolic potential because, beyond the mathematical tractability of the corresponding model, the parabolic profile seems to be suggested by the shape of the actin site, i.e., it appears to be a realist representation of the binding site (for a simplified drawing, see Figure 1). Furthermore, it is also a better choice than that of the piece-wise linear profile, surely less accurate from biophysical point of view.

The stochastic modeling has been confirmed to be one of the more suitable tools for the description and comprehension of such phenomena. Indeed, consider that the Brownian motion $W(t)$ is usually used to represent a thermal agitation (noise) in which such proteins are embedded (see, for instance, ^[7,8,9]). Substantially, as in ^[10], we start by considering a diffusion process for modeling the dynamics of a myosin particle also including time-dependent (forcing) terms in the drift of the diffusion: the $F(t)$ term conveys the action of internal reactions and external loads, causing the tilt of the potential $U(x).$ We call $F(t)$ the tilting force. Such tilt contributes to driving the movement. It is worth mentioning that the inclusion of particular time-dependent tilting forces was born from the need to explain the occurrence of backward Steps of the myosin head (^[10,11]).

Here, our purpose is two-fold: to provide a mathematical generalization of the time-non-homogeneous diffusion model proposed in ^[10] by using the fractional calculus approach for stochastic processes (^[12]), and to provide a more realistic and refined biological model for such phenomenon, in particular, by including memory effects. In such a way, we trust that it will be possible to find a modeling answer to some additional aspects of muscle reactions to causes having a temporal duration, such as, for instance, tiredness, drugs, physical exercise or immobility, or exposition to high or low temperatures. It is clear that such phenomena have a corresponding condition at a microscopic level, i.e., on the actin-myosin interaction that, beyond performing forward-backward Steps, can evolve on a different time scale and embody the above conditions in a correlated behavior on time; in this sense, we say that it has memory. The introduction of time delays can help to improve the modeling of some other biological mechanisms (see, for instance, ^[13,14]).

1.3.   The novelty of the fractional stochastic models

Fractional calculus has proven to be an extremely useful tool for constructing new specialized stochastic processes that are able to model real phenomena including additional physiological evidences (see, for instance, ^[15,16,17]). Here, we show how fractional integrals and derivatives applied for time-changed processes or to stochastic differential equations are able to generate new models with specific peculiarities. The phenomenon considered here, specifically, the possibility that protein interaction takes place on different time-scales and preserves a sort of correlation, can be studied by means of two fractional stochastic modeling processes.

We first propose a time-change process with the long-range dependence property (^[18,19]). The need of applying the time change strategy is justified by the wide approved approach that for many physical phenomena, the time of their evolution proceeds according to the occurrence of other random events (many unknown for now), conditioning the behavior of the phenomenon under observation. Moreover, the adoption of such random time-change induces a dilation of the time and stop periods. The proteins dynamics that appear to not be regular in time and that alternate activity periods and stop periods could be described by this kind of model. Still the actin-myosin dynamics show many obscure sides, and we are confident that such kind of mathematical models can open the door to future understanding such dynamics.

Then, we propose a fractionally integrated process (^[16,20]), whose time-scale can be modulated by the fractional order of a non-local differential operator, involved in its dynamics. Moreover, the motivation of the introduction of a fractionally integrated process is the mathematical nature of this process that is constructed as a solution of a time-fractional differential stochastic equation; in this case, the involved derivative is the integral over time of previous variations of the system we are describing. This means that we are recording the history of the evolution of the system, suitably weighted, and we are using this in the corresponding differential equation. Consequently, the protein's dynamics that embodies all previous variations, more or less long (i.e., with more or less memory) could be described by this kind of model.

1.4.   Dwell times and first passage times

Under additional assumptions, such as that the dynamics is over-damped, meaning that the inertial force is disregarded, we limit to describe the dynamics confined into one space-period of the periodical morphology of the actin filament. Indeed, we assume that the diffusive dynamics occurs in a symmetric one-well parabolic potential, representing the binding site, tilted by a time-dependent force.

Moreover, by taking into account that the potential has a parabolic profile, we model the myosin Steps as it escapes from the potential well: the escape happens when the particle, starting from the bottom of the potential well, crosses over a local maximum of the potential. Note that the escapes can occur in a forward or backward direction. In our models, the processes represent the time-evolution of the myosin head confined in a strip, having two boundaries: the upper with respect to the initial position of the process) and the lower one. A forward Step occurs when the process reaches the upper boundary; a backward Step occurs when the process crosses the lower boundary (see Figure 1).

The dwell time is particularly important in the description of this phenomenon. It represents the time elapsed in the binding site, which, in the modeling context, corresponds to the time elapsing in the potential well before escaping from it. Its investigation can provide useful indications about reaction times and velocity of the muscle contraction.

We model it as the first exit time (FET) of the considered processes from the strip, assuming the boundaries as constant functions in time. In particular, the FET is equal to first passage time (FPT) through the upper boundary or to FPT through the lower boundary. Rigorous definitions of such random variables will be given. Consider also that the attainment of one of two boundaries can be viewed as the variation of the potential value allowing an escape. Moreover, the width of the strip is assumed to be equal to the average distance $L$ between two consecutive binding sites. The initial position of the particle is in the middle of the strip, corresponding to the bottom of the potential well. See Figure 1 as a graphical illustration of such modeling elements (or, for a more detailed illustration, see also Figure 1 of ^[10]).

1.5.   Some previous results in modeling this protein's interaction

In ^[7], by combining loose coupling (^[21,22]) and lever-arm theories (^[4]), we provided a stochastic washboard potential model including biased thermal effects and the power stroke action. We were able to fit the set of data measured by the Yanagida group (^[9]) on the sliding of myosin II heads on immobilized actin filaments for different load conditions. In ^[23], we considered the jump diffusion motion of a particle in a one-well potential; the described dynamics was affected by half-period space shifts at Poisson instants. In ^[24], the case of a double-well potential was also studied. In ^[10], the inclusion of the time-dependent forcing term was considered, and a refined diffusion model was proposed, describing separately the effect of internal and external forces.

1.6.   The main results

In this contribution, by adopting all above modeling assumptions, we start from the diffusion process $X(t)$ with equilibrium level in the bottom of the potential well (in $L/2$), having a time-dependent drift, whose dynamics occurs in the strip $(0, L)$, as proposed in ^[10] (see Section 2). Then, by preserving all involved parameters and functions to represent protein dynamics, we specialize it by applying the stochastic time change by the process $E_\alpha(t)$ (^[25]). The obtained process $X(E_\alpha(t))$ (in Section 3) is a time-changed fractional diffusion process (see, for details, ^[19] and references therein). It is not a Markov process but it preserves the long-range dependence feature of the stochastic time used for the temporal change (^[18]). Furthermore, we specify its mean and covariance functions.

We study its FET probability density function (pdf), particularly useful for understanding of the dwell time (in Section 4). By recalling the subordination operator (see, for instance, ^[26] and references therein), we address the problem by adopting the characterization of the subordinated FPT for the fractional time-changed diffusion processes given in ^[19]. Finally, we are able to provide a new equation for the FET of such processes, specializing, for this case, the Volterra-type integral equation approach (see Proposition 4.1). The results of this study suggest that the time-changed model can be suitable to describe randomly delayed dynamics and also correlated for long times, such as those that can be involved in muscle contraction subject to some illnesses or drugs.

A second fractional model we propose is based on a fractional version of the stochastic differential equation (SDE) of ^[10]. This is obtained by substituting the classical derivative with the fractional Caputo derivative (see, for instance, ^[16,27,28], specialized for neuronal activity models). Indeed, the non-local derivative (see ^[29,30]) allows the construction of a process that preserves the memory of its time evolution. In Section 5, we provide in Proposition 5.12 the explicit expression of the process $\mathcal X(t)$, that is properly a fractionally integrated process (^[20]). Under the specific setting for parameters and involved functions, it is suitable to model the actin-myosin interaction. We also provide its main moments. Some essentials are also recalled for specifying the theoretical setting related to solutions of fraction SDEs. Regarding this model, the study results conclude that it is suitable to describe proteins dynamics strongly sensible to the application of forces, which could correspond to more reactive (sensible) muscle contractions, even if with a greater variability.

In order to prove the applicability of the proposed models, in Section 6, we devised two simulation algorithms: Algorithm 1, for the time-changed process, and Algorithm 2, for the fractionally integrated process. Specifically, Algorithm 1 exploits the transformation formula (3.17) by using the time-changed Brownian motion, whereas Algorithm 2 implements an ad hoc iterative scheme (6.2) based on a discretization formula for the Caputo-fractional derivative. In Section 7, we show and discuss some numerical and simulation results, such as the plot of mean functions, 3-dimensional plot of the covariance, plot of sample paths, and density estimations for simulated FETs. Some concluding remarks complete the manuscript.

2.   A stochastic model for the actin-myosin interaction: the integer case

We start by considering the stochastic model proposed in ^[10] for modeling the actin-myosin interaction, based on the following SDE:

2.1.   Some modeling justifications

The above model born by the corresponding Langevin equation

is widely adopted for describing the over-damped motion of a particle subject to a tilted potential $V(x) = U(x)-Fx$, with $U(x)$ the potential and $F$ the constant tilting force. Specifically, we list all above functions and parameters and the corresponding modelling as follows :

$F,$ the tilting force: it is the sum of an internal force $F_i$ and an external load $F_e$, i.e., we have $F = F_i-F_e$, $U(x),$ the potential assumed to be a space periodic function with period $L$ (as in ^[10], we assume $L = 5.5$ nm equal to the myosin Step size).

$\beta,$ the drag coefficient equal to $90$ pN ns/nm,

T = 293 K, the environmental temperature,

κ_B, the Boltzmann constant.^*

^*Note that $(\prime)$ is used to indicate the derivative of the function with respect to its own argument; in particular, in (2.2), $(\prime)$ is for the space derivative and $(\cdot)$ denotes the time derivative.

In order to include time-varying effects of the tilting force, in ^[10], a time-depending force $F(t)$ was considered in the SDE (2.1) that corresponds to the Langevin equation (2.2) by choosing a parabolic potential $U(x)$ such that

where $U_0$ is the depth of the potential well (see Figure 1 of ^[10]). Hence, due to the form of $V(x) = U(x)-F(t)x,$ we also have

Hence, by setting the parameter $\theta$ as follows

we finally obtain

with W being the standard Brownian motion, and the initial condition is $X(0) = x_0$. For the chosen functions and parameters, the last equation is the same of (2.1): it constitutes the corresponding It$\hat{\text o}$ version of the Langevin equation (2.2). A time-non-homogeneous Ornstein-Uhlenbeck (OU) process $X(t)$ is the solution of SDE (2.1), which exists under suitable assumptions on the tilting force ^[31]. The theory and properties of Gaussian diffusion (GD) (^[32,33]) processes can be exploited for the considered case of equation (2.1).

Note that (2.1) describes the stochastic dynamics of a particle, whose position is described by the value of the stochastic process $\{X(t), t\geq 0\}$ solution of (2.1).

2.2.   The dwell time

By setting the starting point $x_0 = L/2$ (that means in the well of the potential), it is possible to model the dwell time of the myosin as the first exit time (FET) random variable

through the constant boundaries located at the origin (the lower one) and at $L$ (the upper one). The FET pdf is:

The random variable dwell time plays a key role in modeling the sliding of the myosin head along the actin filament and the consequent production of ATP energy. The study of the dwell time and of its distribution is fundamental for the description and comprehension of the mechanism triggering muscle contraction. Indeed, such study is also extremely useful to predict pathological effects of drugs or diseases on muscles.

Beyond a small number of cases for which the pdf of FET is known in closed form, at the moment it is possible to construct approximations of such pdf. For the stochastic process $X(t)$, evaluations of the pdf of the first passage time through a constant level can be obtained by solving with numerical procedures an integral equation (see ^[33]) or by some transformation methods or simulation techniques. In what follows, we give specific details for the considered cases.

3.   A time-changed stochastic model

In order to construct a time-changed fractional model, we adopt the stochastic process $\{X(t), t\geq 0\}$ solution of (2.1) as the parent process, and we compose it with a stochastic time. The time-change is obtained by substituting the time with the positive non-decreasing process $E_\alpha(t)$ that is the inverse of an $\alpha$-stable subordinator process $\sigma_\alpha(t)$. We finally consider in place of $X(t)$ the following process:

By adopting such a model many theoretical results of ^[19] can be exploited and contribute to enrich the description of this biological phenomenon. With this in mind, we recall some mathematical essentials about the time-changed fractional process $Y_\alpha.$

3.1.   Essentials on processes for the time-change

For $\alpha \in (0, 1)$, an $\alpha$-stable subordinator $\sigma_\alpha(t)$ (^[25,34]) is a strictly increasing positive Lévy process that for $\lambda > 0, t > 0$ has the Laplace transform:

with Laplace exponent $\lambda^{\alpha}$ and $\mathbb{E}$ being the expectation operator. A particular property of the subordinator process is the scaling property $\sigma_{\alpha}(t) \stackrel{d} = t^{1/\alpha}\sigma_{\alpha}(1),$ which is extremely useful in simulation algorithms. Note that the notation $\stackrel{d} =$ establishes the equality between finite dimensional distributions of the involved processes.

The inverse $\alpha$-stable subordinator $E_\alpha(t)$ is defined as follows

For any $t > 0,$ the random variables $\sigma_\alpha(t)$ and $E_\alpha(t)$ are absolutely continuous (^[25]). Then, the Laplace-Stieltjes transform of $E_\alpha(t)$ (see, for instance, ^[35,36]) is the following Mittag-Leffler function

where the one-parameter Mittag-Leffler function is defined as

The inverse of $\alpha$-stable subordinator $E_\alpha(t)$ is a self-similar processes, i.e.,

The $E_\alpha(t)$ process has mean

and the covariance function (cf. ^[18] and ^[19]) for $0 < s\leq t$ :

Equation (3.6) involves special functions: the Beta function $B(a, b),$ the hypergeometric function (^[37])

and the incomplete beta function $B(a, b; x)$, which for $x \in [0, 1],$ is defined as

with $B(a, b) = B(a, b; 1).$ It is immediately available from (3.6) the expression of the variance:

It is important to recall the asymptotic power law behavior of the covariance (^[18]), that is

by which the long-range dependence of such a process turns out.

If we denote by $\nu_\alpha(x, t)$ the probability density function (pdf) of $E_\alpha(t)$ and by $\gamma_\alpha(x)$ the pdf of $\sigma_\alpha(1)$, a relationship between them can be highlighted (^[25])

Recalling that $E_\alpha(t)$ assumes positive values for any $t > 0$, its density is zero for $x < 0$, with a discontinuity in $x = 0$ (see the Appendix of ^[19] for further details). In particular, we keep in mind that $E_\alpha$ is an increasing, continuous process with constant values corresponding to the jumps of $\sigma_\alpha.$ Finally, we take note of the Laplace transform of $\nu_\alpha(x, t)$ with respect to $t$:

3.2.   The time-changed process

The time-changed process is the composition of two processes: the parent process $X(t)$ and the inverse of an $\alpha$-stable subordinator $E_\alpha(t)$, independent on $X(t)$. The resulting time-changed process $Y_\alpha(t) = X(E_\alpha(t))$ has continuous sample paths because the parent process is continuous. Specifically, $Y_\alpha(t)$ is constructed by using a time-non-homogeneous GD process as the parent process (the solution process of (2.1)), and the inverse of an $\alpha-$stable subordinator processes for the time-change with pdf $\nu_{\alpha}(s, t)$.

In general, if $f(x, s)$ is the pdf of the parent process, the time-changed process has the following pdf:

We note that $f_{\alpha}(x, t)$ is also called the subordinated density of $f(x, t)$ by means of $\nu_{\alpha}(s, t)$. Theory, properties, and applications of subordinated processes can be found in ^{[29,30,34,38,39,40]}.

Furthermore, by means of the change of variable $ts^{-1/\alpha} = w$ and (3.9), (3.11) can alternatively be written as

In the specific case of the process $X(t)$ solution of (2.1), we know its pdf $f(x, t)$ is a normal density with mean:

and the covariance function

Due to the Gauss-Markov nature of process $X(t)$, its covariance is the product of two functions

whose ratio $r(t) = \rho(t)/\eta(t)$ is a monotonically non-decreasing function.

Hence, from (3.11) or (3.12), we are able to specify the subordinated pdf of the time-changed process $Y_\alpha(t)$. In particular, we can also determine its mean and covariance functions by recalling that for the GD process $X(t),$ the Doob transform holds, i.e.,

where $\{W(t), t\geq 0\}$ is the standard Brownian motion, $r(t) = \rho(t)/\eta(t)$ with $\eta(t), \rho(t)$ being the functions in (3.15) (see ^[33]). From the Doob transform, the time-changed $Y_\alpha$ process is such that

with $\eta(t), \, m_X(t) \in C^1(I)$, $r(t)$ positive monotone increasing $C^1(I)-$function (with $r(0) = 0$) and $W(r(E_\alpha(t)))$ is the time-changed Brownian motion.

Remark 3.1. We remark that it is possible to consider the two time-changed Brownian motions. Indeed, note that the process $W(r(E_\alpha(t)))$ has the following subordinated pdf:

with $f_W(x, r(s))$ the pdf of a standard Brownian motion $W.$ On the other hand, it is also possible to consider another time-changed Brownian motion defined as $W_a(r(t)) = W(E_\alpha(r(t))),$ whose pdf is the following one:

Now, we put into evidence that:

where we applied the substitution of the variable $y = r(s);$ in such a case, the two time-changed Brownian motions are equivalent.

Remark 3.2. (The time-changed Brownian motion $W_\alpha(t) = W(E_\alpha(t))$ and the Caputo derivative)

It is well-known (^[25], ^[41], ^[42]) that the subordinated pdf

where $f_{W}(x, t|y, s),$ the transition pdf of $W(t)$, satisfies the following fractional Fokker-Planck equation

with the initial condition

where $\delta(\cdot)$ is the delta function. In (3.18), the differential operator $D^{\alpha}_{t }$ is the Caputo fractional derivative with respect to $t$, that for a $C^1$ function, $f(t)$ is as in ^[25]:

Due to this fact, the time-changed Brownian motion is also called anomalous or fractional diffusion. Such terminology is also extended to the process $Y_\alpha(t)$, (see ^[19] for details).

3.3.   On the moments of the time-changed process

Note that the time-changed process $Y_\alpha$ is no more Gaussian and no more Markov process.

3.3.1.   The expectation

By exploiting the independence of $E_\alpha$ from $W$, we can calculate the mean of $Y_\alpha$:

From (3.13), (3.5) and (3.2), we finally have:

Such a mean function, in case of a constant force $F(t)\equiv F,$ immediately becomes:

In case of a decaying exponential force $F(t) = Fe^{-t/c}$ with $t, c > 0,$ we have

with $\delta = (1/c-1/\theta),$ from which it finally follows:

It appears evident that the mean is written in terms of several Mittag-Leffler functions.

3.3.2.   The covariance

Then, by recalling the independence of the Brownian motion $W$ from the inverse subordinator $E_a,$ the covariance of $Y_{\alpha}(t)$ can be evaluated, for $s < t,$ as follows:

By taking into account that $cov(W(r(E_\alpha(s))), W(r(E_\alpha(t)))) = \mathbb E(W(r(E_\alpha(s))), W(r(E_\alpha(t)))) = \mathbb E[\min\{r(E_\alpha(s)), r(E_\alpha(t))\}],$ one has

due to the function $r(\cdot)$ and the process $E_\alpha$ both increasing. The explicit expressions of the functions $\eta(\cdot)$ and $r(\cdot)$ of the parent process $X(t)$ are that of the $\eta(t)$ function as in (3.15) and

respectively. Hence, in our specific case, we have

For $t = ks,$ with $k > 0,$ we apply the self-similarity property (3.4) of the process $E_\alpha$, that gives $E_\alpha(t) = E_\alpha(ks) \stackrel{d} = k^{\alpha} E_\alpha(s),$ and finally we obtain

where we used (3.2). The long-range dependence of the covariance of $Y_\alpha(t)$ can be deduced from the asymptotic behaviors of the above Mittag-Leffler functions involved in its expression. Indeed, it is sufficient to recall from ^[43] that, for $0 < \alpha < 1$, and $z > 0,$

This also means that the long-range dependence property of $E_\alpha(t)$ is inherited in the time-changed process $Y_\alpha(t)$.

From (3.32), we can also specify the variance:

4.   The dwell-time pdf by subordinated FET of the time-changed process

In order to clarify the results about the investigation of first passage times of time-changed processes and how these can be used for the specific model of the dwell-time, we resume some essentials. From ^[19], we recall that estimations of the FPT pdf for the time-changed process can be obtained by numerical resolution of an integral equation for the subordinated first passage time through a constant level. Then, we formulate the problem for the specific FET for the time-changed process $Y_\alpha(t)$ and specialize the integral equation approach.

4.1.   About FPT probability density

Indeed, we recall that for a GM process (^[33], ^[10]) with the Doob representation (3.16) and $f[x, t|y, \tau]$ its normal transition pdf, the pdf of the FPT through the level $L$, for $x_0 < L,$ is defined as

with density $g_{T_L}(t|x_0, 0)$. For $m(t), \eta(t), \rho(t), r(t),$ $C^1(I)$-functions defined in subsection 3.2, $g_{T_L}(t|x_0, 0)$ is a solution of a non-singular second-kind Volterra integral equation:

where

For $Y_\alpha(t)$ process, consider the $\alpha$-stable subordinated pdf

Then, denote by $T_{L, \alpha}$ the subordinated FPT with pdf $g_{T_L, \alpha}(t|x_0, 0)$. In Theorem 4.1 of ^[19], was proved that $g_{T_L, \alpha}(L, t|x_0, 0)$ satisfies the equation

with the following integral operator

and

Numerical quadrature procedures applied to equation (4.5) provide approximations of $g_{T_L, \alpha}$ for which no closed form results are available.

4.2.   About the FET pdf

From ^[10] and ^[44], we give some known results about the FET $T_X$ of $X(t)$ defined in (2.6) evolving in the strip $[0, L]\in \mathbb{R}$ with $L > 0$, starting from $X(0) = x_0\in [0, L]$ (in particular, we chose $X_0 = L/2$). It is such that $T_X = \min\{\mathcal T_0, \mathcal T_L\}$, where $\mathcal T_0$ and $\mathcal T_L$ are the FPTs through the lower constant boundary $0$ and the upper one $L$, respectively, i.e.,

Moreover, $g_{\mathcal T_0}(t|x_0, 0)$ and $g_{\mathcal T_L}(t|x_0, 0)$ are the density functions of FPTs $\mathcal T_0$ and $\mathcal T_L$, respectively. Note that the random variable ${\mathcal T_L}$ is different from $T_L$ defined in (4.1), because ${\mathcal T_L}$ is a FPT of the same process through the level $L$ in the presence of an additional absorbing boundary in zero. And, vice versa, ${\mathcal T_0}$ is the FPT of the $X(t)$ process through the level zero in the presence of the additional absorbing boundary in L. Hence, the FET pdf of $T$ of $X(t)$ defined in (2.6) is such that

From ^[44], $g_{\mathcal T_0}(t|x_0, 0)$ and $g_{\mathcal T_L}(t|x_0, 0)$ are solutions of the two following coupled integral equations:

where

with $f_X[x, t|z, \tau]$ being the transition pdf of $X(t)$ and with $S_1 = 0, \, S_2 = L.$

Note that the FET pdf (4.10) of $X(t)$ can be evaluated by applying the standard numerical algorithms for the resolution of this kind of integral equation as (4.11) or a specific algorithm given in ^[44]. Due to the Gaussian distribution of the $X(t)$ process, only the mean and covariance functions of the process are required to apply the numerical procedure; consequently, the functions $\psi_{i}(t|z, \tau)$ ($i = 1, 2$) involved in the integral equations are also completely specified.

4.3.   About the FET of the time-changed process

For the FET of the time-changed process $Y_\alpha(t),$ we can proceed similarly to the case of subordinated FPT. Specifically, we define the subordinated FPT pdfs:

in such a way that the FET is:

Proposition 4.1. The subordinated FPT pdf $g_{\mathcal T_0, \alpha}(t|x_0, 0)$ and $g_{\mathcal T_L, \alpha}(t|x_0, 0)$ are solutions of the two following coupled equations:

and

Proof. The proof follows the proof of Theorem 4.1 of ^[19].

Indeed, we substitute $g_{\mathcal T_0}(t|x_0, 0)$ and $g_{\mathcal T_L}(t|x_0, 0)$ as in (4.11) in the right-hand side of (4.13) and (4.14), respectively, and we have

All involved functions are $L^1$ and the Fubini theorem can be applied to the integral terms at the right-hand side of equations in (4.18), and we can write

and

where we used

Finally, by using (4.17) in the right-hand sides of (4.19) and (4.21), we obtain the coupled equations (4.16) that complete the proof. □

Numerical procedures can be devised in order to solve the coupled equations (4.16), but note that for each value of $t,$ a quadrature procedure is also required for evaluations of (4.17) and finally to obtain approximations of the FET density (4.15). Additionally, evaluations of FPT pdfs for the parent process $X(t)$ are also required. Even if it is a practicable method, a particular effort can be done to reduce the computational cost of such procedures. This will be the object of our future work. An alternative method is provided by simulation algorithms (see Section 6).

5.   Fractionally integrated processes to model the actin-myosin dynamics

Another way to model memory effects in the dynamics of the actin-myosin is to reconsider equation (2.1) and substitute the classical derivative with the fractional Caputo derivative as defined in (3.20), in such a way that the following fractional SDE can be assumed for the model:

Similar fractional models can be found in ^{[16,27,28,45]} and, specialized for neuronal activity models. Here, by following the approach adopted in ^[16], we first recall the theoretical results of ^[46] and specialize them for this kind of equations. Note that the stochastic fractional differential equation (5.1) will be understood in a sense of equation (5) of ^[46] (to avoid collision with the theory developed in ^[47] and ^[48]). Indeed, by comparing with Eq. (2) of ^[46], this fractional SDE can also be written as follows

with $D = \frac{d}{dt}$ the usual derivative. Regarding such a model, we give the same biological interpretation as explained in Introduction and in Section 2, due to the fact that the RHS of Eq.(5.1) is the same of RHS of Eq.(2.1), but with a difference related to the LHS. Indeed, by focusing on the LHS of Eq.(5.2), it is possible to give the following biological justification: differently from the integer case of Eq.(2.1), the solution process $\mathcal X$ of Eq.(5.2) is a process whose time evolution at time $t$ takes into account its whole past evolution until time $t.$ This is essentially due to the presence of the special integro-differential operator on the LHS of Eq.(5.2). In the application context, this process can be adopted for modeling biological phenomenon including some memory (i.e., more or less memory depending of the $\alpha$ value). This could be the case of muscle contraction, and all underlying protein dynamics, in tiredness condition or aging.

5.1.   Some theoretical results about fractional SDE

Equation (5.1) is a Caputo-fractional SDE. Indeed, consider $T > 0$ a real number and the bounded interval $[0, T],$ a typical Caputo-fractional stochastic differential equation is

with $\{W_t\}_{t\in [0, \infty)}$ a real standard Brownian motion in $(\Omega, \mathcal F, \{\mathcal F_t\}_{t\in [0, \infty)}, \mathbb P)$ a complete filtered probability space, $\alpha\in (1/2, 1)$, $A, B, \sigma$ real-valued measurable functions on $[0, T]$ satisfying the following assumptions:

(ⅰ) there exist $L > 0$ such that $\forall x, z\in \mathbb R, \, \forall t\in [0, T]$

(ⅱ) $\int_0^T \Vert B(s, 0)\Vert^2 ds < \infty, \quad$ ess$\sup_{s\in [0, T]} \Vert \sigma(s, 0)\Vert < \infty.$

We are in the suitable space $\mathbb L^2(\Omega, \mathcal F_t, \mathbb P)$, for any $t\in [0, \infty),$ that is the space of mean square integrable functions $f:\Omega \rightarrow \mathbb R$ with the mean square norm $\Vert f \Vert_{ms} = \sqrt{\mathbb E(\Vert f \Vert^2)}.$ Consider a process $\mathcal X: [0, \infty)\rightarrow \mathbb L^2(\Omega, \mathcal F, \mathbb P)$ that is $\mathbb F$-adapted if $X(t)\in \mathbb L^2(\Omega, \mathcal F_t, \mathbb P), \forall t\geq 0.$ Similarly to ^[49], we specify that a classical solution of (5.3) is a $\mathbb F$-adapted process $\mathcal X$ with initial condition $\mathcal X(0) = \mathcal X_0\in L^2(\Omega, \mathcal F_0, \mathbb P)$ such that for $t\in [0, T]:$

or, equivalently,

where $\mathcal I^\alpha$ is the Riemann-Liouville fractional integral (see, for instance, ^[26], ^[50]) of order $\alpha$ defined as follows

with $\Gamma$ the Gamma Euler function, i.e., $\Gamma (z) = \int _0 ^ { + \infty } t^{z-1} e ^{-t} dt \, \ z > 0.$

Moreover, in (5.5), the operator $\mathfrak I^\alpha \left(\sigma(s, \mathcal X(s))dW_s\right)$ is a stochastic fractional integral defined as

We note that this integral can be viewed as a generalization of a fractional integral that defines the Lévy fractional Brownian motion, i.e.,

where $W$ is a standard Brownian motion and the integral has to be interpreted in the Itô sense. We recall that if $\alpha \in \left(\frac{1}{2}, 1\right)$, for fixed $t > 0$, the function $(t-\tau)^{\alpha-1}$ is in $L^2(0, t)$. Thus, the above process is well-defined and it is a Gaussian process starting from 0 at time 0 with probability 1.

Specifically, we recall that in ^[49], it was proved that for any given $X_0,$ there exists a unique solution, whose explicit form was found with a variation constant formula. We can write the specific result in the next proposition.

Proposition 5.1. The solution process of the fractional SDE (5.3) with initial condition $\mathcal X(0) = \mathcal X_0\in L^2(\Omega, \mathcal F_0, \mathbb P)$, $\forall t\in [0, T],$ is:

with $\mathscr E_{\alpha, \alpha}(z) = \sum_{k = 0}^{\infty}\frac{z^k}{\Gamma(\alpha k+\alpha)}$ and $\mathscr E_{\alpha}(z) = \mathscr E_{\alpha, 1}(z)$ Mittag-Leffler functions (^[51]).

The proof is given in ^[46].

5.2.   A fractionally integrated model

We will apply the following corollary to the fractional SDE (5.1)

Corollary 5.1. The fraction SDE:

with $B \in \mathbb L^2([0, T], \mathbb R),$ $\sigma \in \mathbb L^\infty([0, T], \mathbb R)$ admits the following solution process

Note that for $\alpha = 1,$ the fractional derivative and integral reduce to the corresponding classical ones (^[52]). Consequently, the solution process is $X(t)$ of the integer stochastic model.

Proposition 5.2. The solution process of (5.1) is

Proof. By applying Corollary 5.1 to equation (5.1) with the following setting

with all above parameters and functions under assumptions (ⅰ) and (ⅱ) of this subsection, we obtain the solution (5.12).

□

5.3.   Moments of $\mathcal X(t)$

The solution process (5.12) is a Gaussian process (see, for instance, ^[53]:Theorem 4) with the following mean:

The calculus of its covariance leads to:

where the stochastic fractional integral $\mathfrak I^\beta$ is that defined in (5.7). Finally, by solving Eq.(5.14) as in ^[45], we obtain, for $u < t,$

From the covariance we can also derive the variance as follows:

We remark that the provided expressions of the mean (5.13) and of the covariance (5.14) of the process $\mathcal X(t)$ are extremely useful to obtain simulations of these processes. Details of a possible simulation strategy are given in the next section.

6.   Simulation algorithms for fractional stochastic models

We propose some possible simulation algorithms for fractional stochastic models suitably specified for the actin-myosin dynamics. First, we provide the algorithm for the simulation of trajectories of the time-changed process $Y_\alpha(t)$ and the sampling of its first exit times. Then, we propose a further simulation algorithm for the fractionally integrated process $\mathcal X(t)$ of Section 5. Finally, we suggest for comparison another discretization scheme that can be adopted for time-changed fractional stochastic differential equation due to ^[54].

6.1.   Trajectory simulation of the time-changed process

We refer to the time-changed process $Y_\alpha(t)$ studied in Section 3. Trajectories of the time-changed process can be simulated following the Steps listed below, and then samples of simulated dwell-times can be provided. Consequentially, from simulations, it is possible to give estimations of the pdf of the dwell-time of the actin-myosin interaction for the proposed model.

Algorithm 1

The main Steps of the simulation algorithm are:

● Step1: in INPUT, provide the values of parameters: $\alpha, \beta, \theta, L,$ the force $F,$ the Step size value $\Delta t$ for the time discretization, $N$ as the maximum number of time Steps, and $M$ the sample size of dwell times to be simulated;

● Step2: generate in correspondence to time instants $0 = t_0 < t_{\Delta t} < t_{2\Delta t} < \cdots < t_{N\Delta t}$ the random variables $E_\alpha(t_{\Delta t}) < E_\alpha(t_{2\Delta t}) < \cdots < E_\alpha(t_{N\Delta t})$ by using ad hoc functions of R programming packages. In short: use the R-package stabledist to generate the random time $\sigma_\alpha(s)$ of the $\alpha$-stable subordinator; then, according to the definition (3.1), determine the value $E_\alpha(t)$ as the FPT of $\sigma_\alpha(s)$ by $t$.

● Step3: generate a random path of the time-changed Brownian motion in the random times $r(E_\alpha(t_{\Delta t})) < r(E_\alpha(t_{2\Delta t})) < \cdots < r(E_\alpha(t_{N\Delta t}))$, i.e., $W(r(E_\alpha(t_{\Delta t}))), \ldots, W(r(E_\alpha(t_{N\Delta t}))$ as

for $i = 1, \ldots, N$, with $W(r(E_\alpha(t_0)) = 0$, with $Z_i\sim \mathcal N(0, 1)$ as standard normal random variables;

● Step4: evaluate a random path of $Y_\alpha(i\Delta t) = X(E_\alpha(t_{i\Delta t}))$ by means of (3.17) for $i = 1, \ldots, N$, i.e.,

for $i = 1, \ldots, N,$ with function $m_X(\cdot)$ being given in (3.13) and $\eta(\cdot)$ in (3.15);

● Step5: check if $X(E_\alpha(t_{i\Delta t}))$ is over the level $L$ or under level zero; if so, and if it is the first time this occurs, record the correspondent time instant $E_\alpha(t_{i\Delta t})$ that will be an instance of the dwell time (i.e., the FPT of the time-changed process);

● Step6: go to Step2 and repeat the procedure $M$ times with different seeds for the generation of the sequences of $N$ random instants $E_\alpha(t_{\Delta t}) < E_\alpha(t_{2\Delta t}) < \cdots < E_\alpha(t_{N\Delta t})$ providing the $M$ simulated trajectories of the time-changed process.

● Step7: plot an histogram and/or an approximating density for the sample of the simulated $M$ dwell times.

Note that, for the numerical evaluation of Mittag-Leffler functions in Step4, we use the $MittagLeffleR$ R-package.

6.2.   On simulations of the fractionally integrated process

Algorithm for the simulation of fractionally integrated processes $\mathcal X(t)$ of Section 5, we can indicate that it is possible to adopt the simulation algorithm as that in subsection 4.1.1 of ^[20] based on the Cholesky decomposition method of the covariance matrix of the process by which pseudo-Gaussian sample paths can be obtained. Note that in order to apply such algorithm, the numerical evaluation of the mean (5.13) and of the covariance (5.14) of the process is required. High computational costs and possible round-off errors can affect this method.

A second strategy for simulations of sample paths of $\mathcal X(t)$ is that adopted in Section 3 of ^[27], substantially based on the following discretization formula of the Caputo derivative defined in (3.20), for $t = N \Delta t:$

From (6.1), an iterative discretization scheme of the fractional SDE (5.1) can be derived and, specifically, the simulation of sample paths is obtained as follows:

Algorithm 2

● Step1: in INPUT, provide the values of parameters: $\alpha, \beta, \theta, L,$ the force $F,$ the Step size value $\Delta t$ for the time discretization, $N$ the maximum number of time Steps, and $M$ the sample size of dwell times to be simulated;

● Step2: in correspondence to time instants $0 = t_0 < t_{\Delta t} < t_{2\Delta t} < \cdots < t_{N\Delta t}$ construct the path of $\mathcal X(t)$ process adopting the following iterative scheme:

with $Z_{n}\sim \mathcal N(0, \Delta t)$ normal random variables, $\mathcal X(0) = x_0$ and

,

● Step3: check if $\mathcal X(n\Delta t)$ is over the $L$ value or under zero; if so, and if it is the first time this occurs, record the correspondent time instant $n\Delta t$ that will be an instance of the dwell time;

● Step4: go to Step2 and repeat the procedure $M$ times with different seeds for the generation of the sequences of $N$ independent random pseudo-Gaussian numbers $Z_{n}\sim \mathcal N(0, \Delta t), \, n = 1, \ldots, N,$ providing the $M$ simulated trajectories of the process $\mathcal X(t)$.

● Step5: plot an histogram or/and an approximating density for the sample of the simulated $M$ dwell times.

6.3.   Simulations by applying an Euler-type discretization scheme to a fractional time-changed SDE

Now, by referring to some relevant results obtained in ^[47], ^[54], ^[55], ^[56], we consider the time-changed Brownian motion $W_\alpha(t) = W(E_\alpha(t))$ and the process $\mathscr Z_\alpha(t)$ solution of the following equation:

that we propose as an additional model for the protein dynamics. It is considered for comparisons, validations, and further investigations of the proposed models. Note that the process $\mathscr Z_\alpha(t)$ is such that, for $\alpha = 1,$ $\mathscr Z_\alpha(t)\equiv X(t)$ because the (6.4) reduces to (2.1).

Regarding the well-posedness of the considered stochastic equation, we are under all assumptions to guarantee the existence and path-wise uniqueness of the solution (^[54]). Then, we define the solution of the equation (6.4) as the following one

that is the sum of a stochastic integral with respect to $E_\alpha(t)$ and the time-changed Brownian motion $W_\alpha(t).$ The stochastic integral with respect to $E_\alpha(t)$ is well-defined as a stochastic integral with respect to semi-martingales (see, for instance, ^[54]).

Such fractional differential model allows an alternative simulation algorithm based on an Euler-Maruyama discretization scheme (cf. ^[56]). Indeed, from (6.6), we have

where it is possible to follow Step2 of Algorithm 1 to generate random times $E_\alpha(t)$ and to follow Step3 of Algorithm 1 to simulate the values of the time-changed Brownian motion $W_\alpha(t).$

Remark 6.1. A comparison between $\mathscr Z_\alpha(t)$ and the time-changed process $X(E_\alpha(t)):$ the last one can be assumed as the solution of the following equation:

In order to see the similarities and differences between the two processes, it is possible to note that a possible Euler-Maruyama scheme for a simulation algorithm would be the following

This could be an alternative algorithm for simulation of the process $X(E_\alpha(t)),$ even if some consistency conditions should be investigated about the equivalence between the solution process of (6.6) and the process $X(E_\alpha(t))$ constructed in Section 3. All these and other mathematical refinements will be the subject of a next work.

7.   Some numerical results, simulations and biological interpretation

We perform the suggested simulation algorithms in order to construct sample paths of the proposed processes and we derive the statistical density estimations of the simulated dwell times, as described in our models. No quantitative comparisons and/or fit with experimental data is done, but comparisons with the classical integer case are provided by specifying the innovations of such models in terms of the biological modeling of actin-myosin interactions. The following numerical results are given in graphical form in order to show visible validations of the used algorithms and explain the adequacy of the models.

The set of fundamental parameters is done as indicated in Section 2. Other choices are specified in the caption of each figure.

7.1.   Plots of means, simulated paths and FET estimations for $Y_\alpha$ by Algorithm 1

Figure 2 shows the plots of the mean function of $Y_\alpha(t)$ (solid line) and a corresponding simulated path. Specifically, (a)-(b)-(c) refer to mean function (3.25) of $Y_\alpha(t)$ for different $\alpha$ values in the presence of a constant tilting force $F$, whereas (d)-(e)-(f) refer to the mean function (3.27) of $Y_\alpha(t)$ in the presence of a time dependent force $F(t) = Fe^{-t/c}$. The simulated trajectories are performed by means of Algorithm 1. By comparing the classical integer case (a) with fractional cases (b)-(c) (and, respectively, (d) with (e)-(f)), it is possible to see how the choice of the fractional order affects the mean function of the processes and the corresponding behavior of random paths. Periods of constancy in the evolution time of its trajectories appear more evident and frequent as $\alpha$ decreases, in both cases of a constant tilting force $F$ (Figure 2 (a)-(b)-(c)) and also for a time-dependent tilting force (Figure 2 (d)-(e)-(f)). This feature derives from the stretching effect of the time due to the adoption of the time-change regulated by the fractional order $\alpha.$ In terms of biological modeling, such adaptability feature of the proposed model $Y_\alpha$ reveals to be especially useful to describe (and predict) anomalous periods of muscle rest alternating with periods of movement, circumstances that can be verified in the presence of illness conditions. Just in this last case, a delayed muscle dynamics can be observed. Furthermore, periods of immobility, possibly caused by uncontrolled (i.e., random) events, are encoded by the stochastic time $E_\alpha(t).$

Figure 3 shows the plots of the mean function of $Y_\alpha(t)$ subject to a constant tilting force (3.25) on the left side, and to time dependent tilting forces (3.27) on the middle and right side. For different values of the fractional order $\alpha$ used in the time change, the average of a delayed evolution of the corresponding means appears evident in these plots; the more or less evidence of the delay is regulated just by the values of $\alpha.$ In particular, the behavior of the mean function of $Y_\alpha(t)$ is gradually more delayed for decreasing values of $\alpha.$ The biological interpretation analogous to that of Figure 2 can be done, taking into account that in this figure, we provide mean functions corresponding to further values of $\alpha.$

Figure 4 show the plots of approximated FET probability densities of $Y_\alpha$ as estimates of dwell times. It is possible to see how the effect of the fractional order in the time-changed process is relevant also on the profiles of the FET-estimated density for different $\alpha$ values. It is noteworthy the evident heavy tails of these densities, essentially due to the long-range dependence property of the process $Y_\alpha(t).$ In terms of biological modeling, this means that this model is able to describe slow movements represented by means of the non-zero probability of long dwell times (FETs). These circumstances can occur in the presence of a delayed muscle response to stimulus caused by physical damages or mental diseases.

7.2.   Plots of means, simulated paths, and FET estimations for $\mathcal X$ by Algorithm 2

We proceed by presenting our investigations about the process $\mathcal X$ similarly to what was previously carried out for the process $Y_\alpha$, while highlighting differences and similarities.

Figure 5 shows (on the left) the mean functions of the fractional integrated process $\mathcal X(t)$, their tails (on the middle), and some corresponding simulated paths (on the right) by means of Algorithm 2 for specified $\alpha$ values when a time-depending tilting force is applied. The main feature we captured by this numerical investigations is that such a process reacts to the tilting force; alternatively, we can say that it is very sensitive to the tilting force due to the high value of the mean in correspondence to short times. In particular, the peak of the mean plot is higher for higher values of $\alpha,$ until $\alpha = 0.9$, which is closer to the classical integer case. Instead, successively, i.e., for long times, we can observe a rapid decay differently from the mean of $Y_\alpha$ in Figure 3. Look then at the behavior of the means in the tails in the middle of Figure 5 and observe that an opposite order of tails occurs: indeed, the decay is slower for lower values of $\alpha.$ In a biological context, this can be interpreted as an instantaneous reaction of the mean behavior of the particle to the effect of the tilting driving force similarly to the mean of $Y_\alpha$ (compare with the right side of Figure 3), even if for $Y_\alpha$ such a reaction persists for longer. Such a reaction appears to be less evident for lower values of $\alpha$. We can still specify that, due to the mathematical construction of such a process, the order $\alpha$ of the fractional derivative affects the process by including memory effects in its time evolution as $\alpha$ decreases. This is also shown in the middle side of Figure 5, regarding the behavior of the mean tails. Hence, for lower values of $\alpha,$ the model seems less sensitive to the effect of the driving force, but the model still encodes such effect in longer times, which can be explained as the preservation of memory effects.

Figure 6 shows the estimations of dwell time (FET) density from $10^3$ simulated paths by means of Algorithm 2 of the fractionally integrated process $\mathcal X(t)$. In this case, (quasi) heavy tails also appear in the FET density estimations, but they appear very early, i.e., they are visible on a reduced time period with respect to that of $Y_\alpha$ case. Indeed, in this case, we can talk about a preservation of memory effects as $\alpha$ decreases, while for $Y_\alpha$ we can talk about an enlarged time-scale. We understand that this phenomenon is also due to the increased variability (variance) that is in this model as $\alpha$ decreases. In this regard, compare the simulated paths of $Y_\alpha$ (Figure 2) and those of $\mathcal X$ (in the right side of Figure 5 right): we will show that the difference is also due to the different variance values of the two processes (see subsection 7.3).

7.3.   Numerical evaluations of covariances and comparisons

We investigate how the behavior of the covariance and variance of the two proposed models changes for different values of $\alpha.$

In Figure 7, we plot the covariance of $Y_\alpha(t)$ as in Eq. (3.32) on the left and of the variance as in Eq. (3.33). The resulting curves are ordered by values of $\alpha.$ The covariance curves decay slower for lower values of $\alpha$; we expect this due to the power decay already discussed about this covariance. The variance curves increase slower for lower values of $\alpha.$

A wide investigation into covariance is performed by using 3-dimensional plots; indeed, Figure 8 is composed by three different 3-dimensional plots of covariance for 3 different values of $\alpha.$ Figure 9 completes the illustrations of the $cov(Y_\alpha(s), Y_\alpha(t))$ as in Eq. (3.32) in order to improve the understanding of the behavior of such a function. In particular, the color maps put in evidence that the correlation persists even for long periods especially for low values of $\alpha.$ This feature allows to adopt such a model when the proteins interaction seems randomly delayed (as well as the muscle contraction).

Analogous investigations were carried out for process $\mathcal X.$ Figures 10, 11, 12 correspond to Figures 7, 8, 9. The behaviors of the covariance and variance of $\mathcal X$ are different than those of $Y_\alpha.$ The covariance decay occurs more rapidly, whereas the variance increases as $\alpha$ increases. Then, we can focus on the color maps of Figure 12 in which it is visible that an extreme concentration of high values of the covariance is in the variance and just for $\alpha = 0.6.$ Compare Figure 12 with Figure 9: the high covariance values are distributed in different zones in the plane $(s, t).$

For biological applications, we understand that this feature of the two models allows to distinguish when a model is preferable to the other one. Indeed, if we need to describe the phenomenon enlarging time, we can adopt the $Y_\alpha$ process that preserves the correlated behavior. On the other hand, if we want to preserve the history of the phenomenon, we can use the $\mathcal X$ process that includes memory effects. In this last case, the associated higher variance can sometimes also be used to explain some anomalous events.

7.4.   Final comparisons between $Y_\alpha$ and $\mathcal X$ models

Finally, in Figure 13, we compare the FET density estimations of the two processes $Y_\alpha$ and $\mathcal X$ subject to the same tilting force. Differences: the height of the peaks is greater for $\mathcal X$ (on the right side) than for $Y_\alpha$ (on the left side); hence, the process $\mathcal X$ seems to be more sensitive to the effect of the tilting force in such a way that there is also much more probability that the exits occur quite quickly, before those of $Y_\alpha$. Similarities: in both cases, we note heavier tails as $\alpha$ decreases but longer for $Y_\alpha$. Consequently, one can conclude that the heavy-tailed FET densities come out for both processes, but, after the detailed analysis of the covariances and variances of the two processes, we can conclude that, for the time-changed process $Y_\alpha$, this is a consequence of the long-range dependence property, while in case of $\mathcal X$ process it is related to the increasing variance as $\alpha$ decreases.

The heavy-tailed distributions founded in both models, even if triggered by different causes, turn out especially useful for modeling the actin-myosin dynamics because the corresponding models can explain and predict exits (jumps) occurring after long times and also include (under suitable tilting forces) specified percentages of backward jumps. In order to confirm this last deduction, a quantitative analysis of both models by means of validation procedures on real data will be the object of a future work.

Furthermore, to validate the capability of the proposed model for modeling not only dwell times related to forward Steps of the myosin head but also to backward ones, we show in Figure 14 that for a suitable choice of the values of parameters and the tilting forces (as indicated in the caption of Figure 14) the trajectories with backward Steps (note that these sample paths also attain negative values) can also be generated for time-changed processes (left side of Figure 14) and for fractionally integrated processes (right side of Figure 14).

7.5.   Some concluding remarks

We propose two fractional stochastic models for the description, understanding, and prediction of the mechanism triggering muscle contraction: the actin-myosin dynamics and the corresponding dwell time. We started with a stochastic model ^[10] that was able to fit many experimental features and data, and we generalized such a model by means of the fractional calculus approach. The main motivation of such proposed models was born from the need to provide more refined models including some memory effects.

The first proposed model is based on a time-changed diffusion process, for which we provide some theoretical details. In particular, we are able to adapt the integral equation approach for the corresponding FET from a strip. This model is characterized by a long-range dependent covariance, and such property is due to the applied time-change that stretches, in some sense, the time.

The second fractional model is based on a fractionally integrated process, which is the solution of a Caputo-fractional stochastic differential equation for which we provide expressions for the mean and covariance function. This model has memory of its time-evolution, a consequence of the use of the fractional non-local derivative.

In order to obtain statistical estimations of the dwell time, we provide two algorithmic schemes useful to simulate these kinds of processes. Only to validate the simulation algorithms and the agreement with the theoretical settings, we perform simulations of paths of the considered processes and show how to derive statistical estimations of the corresponding dwell time pdf.

A detailed analysis of the covariances of the two proposed models is provided by means of several numerical evaluations, figures, comparisons, and comments related to the potential applications to the biological context.

For application purposes, we can conclude that the time-changed $Y_\alpha$ model can describe biological randomly delayed dynamics, but correlated for long times, whereas the $\mathcal X$ model proves to be suitable for more impulsive reactions of muscle contractions, associated with more variability.

An ad hoc implementation of such algorithmic schemes for a detailed investigations about forward and backward Steps of the myosin head will be successively carried out. In addition, the adaptability of the proposed models and the simulation approach to real experimental data is the main purpose of our next investigation.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was partially funded by the projects: PRIN-PNRR P2022XSF5H of European Union – NextGenerationEU–(MUR) Project Title “Stochastic Models in Biomathematics and Applications”; PRIN2022-MUR 2022XZSAFN; Project DAISY: “Dynamic Analysis of Interacting Biological Systems through Statistics” of the Univ. Vanvitelli (D.R. 111, 09/02/2024). Pirozzi is a member of GNCSINdaM group.

The authors thank the anonymous reviewers for their precious suggestions and greatly appreciates comments, which led to the improved revised version of this paper.

Conflict of interest

The authors declare there is no conflict of interest.