Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Fundam Res . 2025 Jan 11;6(2):659–671. doi: 10.1016/j.fmre.2024.12.024 Search in PMC Search in PubMed View in NLM Catalog Add to search Multiscale modeling and simulation for anomalous and nonergodic dynamics: From statistics to mathematics Heng Wang Heng Wang a School of Mathematics and Statistics, State Key Laboratory of Natural Product Chemistry, Lanzhou University, Lanzhou 730000, China Find articles by Heng Wang a , Xuhao Li Xuhao Li b School of Mathematical Sciences, Anhui University, Hefei 230601, China Find articles by Xuhao Li b , Lijing Zhao Lijing Zhao c School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an 710129, China Find articles by Lijing Zhao c , Weihua Deng Weihua Deng a School of Mathematics and Statistics, State Key Laboratory of Natural Product Chemistry, Lanzhou University, Lanzhou 730000, China Find articles by Weihua Deng a, ⁎ Author information Article notes Copyright and License information a School of Mathematics and Statistics, State Key Laboratory of Natural Product Chemistry, Lanzhou University, Lanzhou 730000, China b School of Mathematical Sciences, Anhui University, Hefei 230601, China c School of Mathematics and Statistics, Northwestern Polytechnical University, Xi’an 710129, China ⁎ Corresponding author. [email protected] Received 2024 Sep 23; Revised 2024 Dec 24; Accepted 2024 Dec 25; Collection date 2026 Mar. © 2025 The Authors. Publishing Services by Elsevier B.V. on behalf of KeAi Communications Co. Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). PMC Copyright notice PMCID: PMC13069665 PMID: 41971797 Abstract In recent decades, anomalous and nonergodic diffusions are topical issues in almost all disciplines. In 2004, the phrase “anomalous is normal” was used in a title of a Physical Review Letters (PRL) paper, which reveals that the diffusion of classical particles on a solid surface has rich anomalous behavior controlled by the friction coefficient, meaning that anomalous diffusion phenomena are ubiquitous in the natural world. This review article first builds the microscopic models (stochastic processes) to describe the experimentally observed phenomena of the motion of the physical/abstract particles under the frameworks of continuous time random walks, Langevin picture, subordinated strong Markov process, etc. Beyond directly conducting statistical analyses on these stochastic processes, we target on delving into these microscopic models to uncover their physical mechanisms and digging out their potential applications. According to the application scenarios and the research requirements, we design the appropriate statistical observables, e.g., the position of the particles, functional of the trajectories of the particles, probability of the first passage time, escape probability, etc. Then, we derive the governing equations of the probability density functions of the statistical observables. We do the mathematical studies on the equations, including well-posedness and regularity analyses, designing numerical schemes, performing numerical analyses, etc. Finally, we present the applications of the models in chemistry and biology, and propose future prospects in this research field. Keywords: Multiscale modeling, Anomalous and nonergodic dynamics, Statistical observables, PDE analyses, Scientific computation, Deep learning 1. Introduction In nature, “motion” is constantly occurring in a variety of forms. As the process of gene transcription initiates, RNA polymerase moves forward along the DNA synthesizing mRNA and exhibits three different states of motion: transcription extension, backtracking and backtracking recovery [ [1] , [2] , [3] ]. Molecular motors operate akin to efficient couriers, tasked with the precise delivery of the newly synthesized mRNA to target sites [ 4 ]. The movement of neurotransmitter receptors on the postsynaptic membrane displays an anomalous diffusion of the two states as a consequence of the confinement of nanoclusters [ 5 ]. Ligand-receptor interactions have been observed to trigger skipping in three-dimensional media and sliding mechanisms on two-dimensional living cell membranes [ 6 ]. Telomeres, special structures at the ends of chromosomes, act as sentinels guarding against genomic instability, yet they are subject to gradual shortening with each cell division [ [7] , [8] ]. In the processes of polymerization and depolymerization, polymers exhibit Brownian non-Gaussian kinetic characteristics [ 9 ]. In environments where molecules are densely packed, tracer polymers exhibit anomalous diffusion [ 10 ]. The irregular connectivity of pore spaces gives rise to anomalous behavior in fluid flow and chemical transport processes [ 11 ]. In the broad scope of ecology, the movement of animals across extensive spatial and temporal scales frequently exhibits anomalous dynamical characteristics [ [12] , [13] ]. To elucidate the essence of these natural phenomena and uncover the underlying physical mechanisms, scientists employ statistical and mathematical tools to quantify the dynamical behaviors. Methods for quantitative modeling across multiple scales are primarily categorized into two types: one focused on simulating dynamics at the microscale, and the other dedicated to deriving or establishing evolutionary equations at the macroscale. There are two commonly used modeling frameworks, the continuous time random walk (CTRW) model and the Langevin equation. The CTRW model [ 14 ] assumes that the process of particle motion is a wait-jump periodic process, and that the waiting time τ i and the jump length ξ i are two random variables satisfying the joint probability density function (PDF) ϕ ( x , t ) = E [ δ ( ξ i − x ) δ ( τ i − t ) ] . Therefore, the probability densities for the waiting time and the jump length are ψ ( t ) = ∫ − ∞ + ∞ ϕ ( x , t ) d x and η ( x ) = ∫ 0 + ∞ ϕ ( x , t ) d t , respectively. When the waiting time is independent of the jump length, one can obtain the decoupled form ϕ ( x , t ) = η ( x ) ψ ( t ) . For a given ϕ ( x , t ) , the particle trajectory can be plotted as in Fig. 1 . The CTRW model has proven to be a robust and effective tool for quantifying anomalous contaminant transport in porous and faulted geological formations [ 15 ]. The CTRW model with a “long memory” waiting time distribution has demonstrated its effectiveness in the financial sector. Assuming that the reaction waiting times for promoter transitions, mRNA synthesis and degradation follow a general distribution, J. Zhang et al. coarsely derive an iterative equation for the moment-generating function of copy number of mRNA, which conveniently calculates mRNA raw and binomial moments of any order [ 16 ]. E. Roldán et al. examine the distribution of recovery time by describing the kinetic behaviour of RNA polymerase during the retrospective recovery phase as a CTRW [ 17 ]. Subsequently, W. H. Deng et al. expand upon the CTRW model by introducing the multi-internal states modeling approach [ [18] , [19] ] and the alternating states modeling method [ [20] , [21] ], and then also establish the Lévy walk model in non-static media [ 22 ]. S. Fedotov et al. establish a two-state nonlinear CTRW model to depict tumor-cell migration and proliferation invasion [ 23 ]. An intriguing addition to the models is the stochastic resetting process, which mimics the behavior of reverting to the starting point to initiate a new search after an unsuccessful attempt [ [24] , [25] ]. Fig. 1. Open in a new tab Particle trajectory of the CTRW model . Langevin’s equation is an application of Newton’s second law “ F = m a ” to Brownian particles, which has the form [ 26 ] m d 2 x ( t ) d t 2 = − m γ d x ( t ) d t + ξ 1 ( t ) , (1) where x ( t ) is the position of the particle at time t , m is the mass of the particle, − m γ d x ( t ) d t is the frictional force, and ξ 1 ( t ) is the random force (Gaussian white noise, due to the random collisions of the surrounding molecules). When the random force is no longer white, the general form of the Langevin equation can be written as m d 2 x ( t ) d t 2 = − m ∫ t 0 t γ ( t − t ′ ) d x ( t ′ ) d t d t ′ + ξ 2 ( t ) , (2) and the fluctuation-dissipation theorem gives the relationship between the friction coefficient and the random force [ 27 ]. A. D. Viñales et al. introduce a Mittag-Leffler correlated random force leading to anomalous diffusion [ [28] , [29] ]. T. Sandev et al. study the analytical form of the relaxation functions for the three-parameter generalized Langevin equation [ 30 ]. The general form of the Langevin equation, where the friction term is represented by the regularized Prabhakar derivative, is discussed in [ 31 ], and the form incorporating a tempered memory kernel is considered in [ 32 ]. Analyzing the “social forces” felt by pedestrians allows us to model their dynamic behavior as a set of nonlinear coupled Langevin equations [ 33 ]. Reference [ 34 ] presents the Langevin equation for polymers engaged in polymerization/depolymerization reactions, by establishing a random diffusion coefficient that correlates with particle size. Reference [ 35 ] presents the Langevin equations for continuous time Lévy walks. X. D. Wang et al. present the Langevin description of the Lévy walk [ 36 ]. Y. Chen et al. then provide the Langevin description of the Lévy walk with memory [ 37 ] and examine the impact of an external force [ [38] , [39] ]. By making a stochastic scale variation of time, one can construct the time-changed (subordinated) Langevin equation [ [40] , [41] ]. The path properties of the subordinated Brownian motion have been investigated in [ 42 ]. Additionally, Ref. [ 43 ] has considered the time-changed fractional Brownian motion, discussing the moments, correlation structure, and so on. References [ [44] , [45] , [46] ] study the Langevin equations for the subordinated Brownian noise of the tempered Mittag-Leffler memory kernel. Let S α ( t ) denote the strictly increasing α -stable Lévy process (subordinator) with the well-known formula for its Laplace transform 〈 e − k S α ( t ) 〉 = e − t k α . Define E α ( t ) = inf { τ : S α ( τ ) > t } , which denotes the left inverse of S α ( t ) . Then there are scale limit relationships between the CTRW models and the subordinated Brownian motions [ [40] , [47] ], as shown in Table 1 . Table 1. The scale limit relationships between the CTRW models and the subordinated Brownian motions . CTRW model subordinated Brownian motion 〈 τ i 〉 < ∞ , 〈 ξ i 2 〉 < ∞ B ( t ) 〈 τ i 〉 < ∞ , ξ i ∼ | x | − 2 β − 1 B ( S β ( t ) ) τ i ∼ t − α − 1 , 〈 ξ i 2 〉 < ∞ B ( E α ( t ) ) τ i ∼ t − α − 1 , ξ i ∼ | x | − 2 β − 1 B ( S β ( E α ( t ) ) ) Open in a new tab Based on the microscopic model, partial differential equations (PDEs) governing the probability distribution of particle positions can be derived by integral transformations [ [40] , [48] , [49] , [50] ]. It is also possible to derive the PDEs from the macroscopic scale according to physical assumptions. Furthermore, it is also possible to derive equations that describe the probability distributions governing the macroscopic statistical properties of physical processes, which include the position functional of the particle, the probability density of the mean first passage time, and even the escape probability [ [34] , [51] , [52] , [53] ]. For the PDEs mentioned above, researchers also discuss the well-posedness and regularity of the solution, as well as develop the corresponding computational methods. The traditional numerical computational methods, such as the finite difference method, the convolution integral method, and the variational method [ [54] , [55] , [56] , [57] , [58] , [59] ], provide effective and highly accurate approximations for the solutions of PDEs, but their computational complexity grows exponentially with the dimensionality [ 59 ]. With the growth of data resources and computational power, deep learning has become an important tool for solving high-dimensional problems in scientific research. In recent years, a large number of deep learning-based PDE solvers have been developed, most of which are inspired by traditional methods. In 2018, W. E et al. develop the deep Ritz method, which uses deep learning to solve variational problems corresponding to PDEs [ 60 ]. Least squares-based deep learning methods include the deep Galerkin method [ 61 ] and physics-informed neural networks [ 62 ], which train models by minimizing the squared residuals of PDEs. Weak adversarial networks [ 63 ] offer an approach to tackle the weak formulations of high-dimensional PDEs utilizing adversarial learning techniques. J. Han et al. propose a deep learning method for solving parabolic PDEs based on backward stochastic differential equations (BSDEs), called the deep BSDE method [ [64] , [65] ]. References [ [66] , [67] ] provide posterior estimates of the deep BSDE method. Subsequently, H. Wang et al. develop the deep BSDE method for solving infinite-dimensional coupled polymer diffusion systems [ 68 ]. Based on the aforementioned discussions, this paper provides a comprehensive review of the research progress in multiscale modeling of anomalous non-ergodic dynamical systems found in nature and offers insights into future research avenues. The review is structured as follows: In the second section, we delve into microscopic models, derive the PDEs that govern the probability distributions of various statistical observables, and present some theoretical statistical analysis results. The third section is dedicated to the theoretical analysis and computational methods for the macroscopic equations. In the fourth part, we showcase some applications of the developed models in physics, biology, and engineering. Lastly, we propose potential directions for future research prospects. 2. Mathematical modeling This research focuses on the motion of real physical particle or abstract “particle”, for example, the fluctuation of the financial market. So the first step is to use mathematical language to describe the motion, that is, build the microscopic model (stochastic process) to characterize the motion. Three types of microscopic models will be discussed here, i.e., CTRW, Langevin equation, miscellaneous ones, including subordinated Brownian motion, alternating process, and resetting process, etc. Beyond conducting direct statistical analyses on these stochastic processes, we target on delving into these microscopic models to uncover their physical mechanisms, and digging out their potential applications. To accomplish this, we initiate by designing the corresponding statistical observables and deriving the governing equations that dictate the probability distributions of these statistical observables. 2.1. Microscopic models 2.1.1. Continuous time random walks In mathematics, CTRW is a stochastic process with arbitrary given distributions of jump lengths and waiting times, originally discussed by E. W. Montroll and G. H. Weiss [ [14] , [69] ]; and the CTRW is first applied to physical systems by H. Scher and M. Lax [ 70 ], where the random variables of waiting times and displacements are independent and identically distributed, respectively. For the past several decades, CTRW model is used to describe different kinds of anomalous systems, ranging from amorphous semiconductors to DNA molecules. The CTRW model with multiple internal states is built in [ 19 ]. We introduce multiple internal states to characterize some natural phenomena, such as, traps in amorphous semiconductors and ionic currents in cell membranes. Each internal state corresponds to a distribution of waiting time and jump length, and the transitions between internal states are described by a Markov chain with a transition matrix M . The dimension of matrix M is n × n , where n is the number of internal states. The element m i j of matrix M represents the transition probability from state i to j . Lévy walk is a CTRW model with the jump length determined by the generated waiting time, i.e., for the waiting time τ i , the jump length equals to v τ i , in which v is the speed of the particles. The Lévy walk with multiple internal states is built in [ 18 ], having the same structure as the CTRW one. 2.1.2. Langevin equations In physics, “external potential” refers to the potential energy imposed on a physical system by external factors or environment, generally denoted as V ( x ) . The Langevin picture for depicting external potentials offers significant advantages compared to the CTRW model. The dynamical model of a Brownian particle subject to an external potential can be represented in the Langevin equation, that is, m d 2 x ( t ) d t 2 = − ∇ V ( x ( t ) , t ) − m γ d x ( t ) d t + ξ 1 ( t ) , (3) where V ( x , t ) is the external potential and γ is the friction coefficient. If considering the delayed effect of friction, one can use a natural generalization of the Langevin equation with external potential, known as the generalized Langevin equation m d 2 x ( t ) d t 2 = − ∇ V ( x ( t ) , t ) − m ∫ 0 t γ ( t − t ′ ) d x ( t ′ ) d t d t ′ + ξ 2 ( t ) . (4) It is important to note that ξ 1 ( t ) and ξ 2 ( t ) in Eqs. 3 and 4 are different. Since both of them are internal noises, according to the fluctuation-dissipation theorem [ 27 ] (friction and random driving forces come from the same source), their correlation functions are δ ( t ) = 1 m k B T E [ ξ 1 ( t ′ ) ξ 1 ( t ′ + t ) ] and ( ∇ V ( x , t ) ≡ 0 ) γ ( t ) = 1 m k B T E [ ξ 2 ( t ′ ) ξ 2 ( t ′ + t ) ] , respectively, where k B is the Boltzmann constant and T is the absolute temperature. It is observed in [ 71 ] that when colloidal beads diffuse along linear phospholipid bilayer tubes or through entangled F-actinnetworks, the motion of beads show Brownian yet non-Gaussian dynamics. Later, the diffusing diffusivity microscopic model describing this kind of phenomena is built as [ 72 ] d x ( t ) d t = 2 D ( t ) ξ ( t ) , (5) D ( t ) = y 2 ( t ) , (6) d y ( t ) d t = − y ( t ) + η ( t ) , (7) where, x ( t ) and y ( t ) (Ornstein-Uhlenbeck process) are stochastic processes driven by independent Gaussian white noises ξ ( t ) and η ( t ) , respectively, and D ( t ) = y 2 ( t ) provides the diffusion coefficient for x ( t ) . The reasons for choosing the square of the Ornstein-Uhlenbeck process as the diffusing diffusivity of x ( t ) are threefold. Firstly, the non-negativity of this choice avoids the need to add a reflecting boundary condition when D ( t ) = 0 , making the analysis easier to handle. Secondly, it ensures that the dynamics of the diffusivity is stationary under the given correlation time. Moreover, it also guarantees that the PDF of D ( t ) has exponential tails, resulting in a Laplace-like distribution for x ( t ) at short times. Over long times, the particle’s motion with an effective diffusion coefficient 〈 D 〉 s t = lim t → ∞ 〈 y 2 ( t ) 〉 exhibits a crossover to a normal distribution. The Brownian non-Gaussian diffusion induced by polymerization is discussed in [ [9] , [68] ], being modelled by Eq. 5 with the diffusing diffusivity D ( N ( t ) ) , where N ( t ) is the birth-death process satisfying P ( N ( t + τ ) − N ( t ) = k | N ( t ) = n ) = { α ( n ) τ + o ( τ ) , k = 1 , β ( n ) τ + o ( τ ) , k = − 1 , 1 − ( α ( n ) + β ( n ) ) τ + o ( τ ) , k = 0 , o ( τ ) , otherwise , (8) with α ( n ) , β ( n ) ≥ 0 for n ∈ N , and β ( 0 ) = 0 . 2.1.3. Miscellaneous processes In probability theory, the technique of time-changing for Brownian motion is one of the important methods to build a new stochastic process. First, one needs to define a time change process, e.g., α -stable Lévy process S α ( t ) with α ∈ ( 0 , 1 ) , and its inverse E α ( t ) . Then, B ( S α ( t ) ) and B ( E α ( t ) ) , which are the time change to Brownian motion B ( t ) , can be respectively used to describe superdiffusion, being the 2 α -stable Lévy process, and subdiffusion. In fact, one can also build the stochastic process B ( S α ( E β ( t ) ) ) ( β ∈ ( 0 , 1 ) ) to characterize the competition between superdiffusion and subdiffusion. One can build multiple-state alternating stochastic process, for example, a two-state alternating stochastic process [ 21 ], which alternates between Lévy walk and Brownian motion, i.e., the sojourn times in the Lévy walk and Brownian motion stages follow power-law distributions with exponents α + and α − ( 0 < α ± < 2 ) , respectively. This kind of processes generally have several different scales and show strong anomalous diffusion phenomena [ 21 ]. Resetting process is needed in some practical applications, e.g., search processes [ [73] , [74] , [75] , [76] ], stochastic algorithms [ [77] , [78] ], catastrophic phenomena in population dynamics [ [79] , [80] , [81] ], etc. A concrete example is given as follows [ 25 ]: the position x ( t ) of the particle with an initial position of x ( 0 ) = x 0 at time t is updated by the stochastic rule x ( t + d t ) = { X r , with probability r d t , x ( t ) + d B ( t ) , with probability ( 1 − r d t ) , (9) where X r is a fixed position, r is the resetting rate, and B ( t ) is Brownian motion. This stochastic process describes a Brownian motion particle with Poissonian resetting, which, after a time interval d t , is reset to position X r with a probability of r d t , and undergoes Brownian motion with a probability of ( 1 − r d t ) . The motion in non-static media is an important research field. In biology, during the development of vertebrate embryos, cell migrations occur on an underlying tissue domain in response to some factor, such as nutrient. Over the time scale of days in which this cell migration occurs, the underlying tissue is itself growing. In this case, the media is an expanding one with scale factor a ( t ) . Following the discussions from (1)-(3) to (4) in [ 82 ], the original Langevin equation is modified as d x ( t ) d t = d a ( t ) d t x ( t ) a ( t ) + ξ ( t ) . (10) 2.2. Macroscopic equations This subsection presents the governing equations of the probability distributions of some typical statistical observations for the above built microscopic models. 2.2.1. Fokker-Planck equations Here we focus on the probability distribution of the position of the particles [ [40] , [47] ]. The most well-known case is Brownian motion B ( t ) , its Fokker-Planck equation behaves as { ∂ W ( x , t ) ∂ t = Δ W ( x , t ) , W ( x , t = 0 ) = δ ( x ) . (11) If the particle performs subdiffusion, microscopically modelled by B ( E α ( t ) ) or the corresponding CTRW or Langevin picture, the governing equation of the probability distribution of its position is described as [ 47 ] { ∂ W ( x , t ) ∂ t = 0 D t 1 − α Δ W ( x , t ) , W ( x , t = 0 ) = δ ( x ) , (12) where 0 D t 1 − α is the Riemann-Liouville operator defined as 0 D t 1 − α W ( x , t ) = 1 Γ ( α ) ∂ ∂ t ∫ 0 t W ( x , t ′ ) ( t − t ′ ) 1 − α d t ′ . Similarly, the superdiffusion of the particle can be microscopically modelled by B ( S β ( t ) ) , or through the corresponding CTRW or Langevin picture. Its governing equation of probability distribution of the position is given as [ [83] , [84] ] { ∂ W ( x , t ) ∂ t = − ( − Δ ) β W ( x , t ) , W ( x , t = 0 ) = δ ( x ) , (13) where F { − ( − Δ ) β W ( x , t ) } ( k ) = − | k | 2 β F { W ( x , t ) } ( k ) . More general case including possible normal diffusion, subdiffusion, and superdiffusion can be microscopically described by the competition model B ( S β ( E α ( t ) ) ) . The corresponding macroscopic equation behaves as [ [85] , [86] ] { ∂ W ( x , t ) ∂ t = 0 D t 1 − α ( − ( − Δ ) β ) W ( x , t ) , W ( x , t = 0 ) = δ ( x ) . (14) There are also more generalized models of the space-time fractional diffusion equations by introducing general memory kernel [ [87] , [88] , [89] ]. The forward and backward Fokker-Planck equations for the polymer diffusion with the microscopic model Eq. 5 and Eq. 8 are respectively described as [ [34] , [68] ] { ∂ W ( n , x , t ) ∂ t = L n W ( n , x , t ) + D ( n ) ∇ x 2 W ( n , x , t ) , W ( n , x , t = 0 ) = g ( n ) δ ( x ) , (15) and { ∂ W n 0 ( x , t ) ∂ t = F n 0 W n 0 ( x , t ) + D ( n 0 ) ∇ x 2 W n 0 ( x , t ) , W n 0 ( x , t = 0 ) = δ ( x ) , (16) where L n and F n 0 are given difference operators defined as L n W ( n , x , t ) = { μ ( n + 1 ) W ( n + 1 , x , t ) − μ ( n ) W ( n , x , t ) + λ ( n − 1 ) W ( n − 1 , x , t ) − λ ( n ) W ( n , x , t ) , n > 0 , μ ( 1 ) W ( 1 , x , t ) − λ ( 0 ) W ( 0 , x , t ) , n = 0 , and F n 0 W n 0 ( x , t ) = { λ ( n 0 ) ( W n 0 + 1 ( x , t ) − W n 0 ( x , t ) ) + μ ( n 0 ) ( W n 0 − 1 ( x , t ) − W n 0 ( x , t ) ) , n 0 > 0 , λ ( 0 ) ( W 1 ( x , t ) − W 0 ( x , t ) ) , n 0 = 0 . 2.2.2. Feynman-Kac equations We turn to the statistical observable A , which is the functional of the trajectory of a particle x ( t ) , defined as A = ∫ 0 t U ( x ( t ′ ) ) d t ′ , with U ( x ) being a given function. Two types of Feynman-Kac equations have been derived, respectively governing the joint distribution of A and x ( t ) (called forward Feynman-Kac equation) and the distribution of A with the process’s starting point x ( 0 ) = x 0 as the parameter (called backward Feynman-Kac equation). First, we present the backward Feynman-Kac equation for the microscopic model with external potential and multiplicative noise. The model behaves as d x ( t ) d t = − ∇ V ( x ( t ) ) + g ( x ( t ) ) ξ ( t ) . (17) While ξ ( t ) can be white noise and β -stable Lévy noise ( g ( x ) ≡ 1 ), the backward Feynman-Kac equations of Eq. 17 are respectively given as [ 52 ] { ∂ G x 0 ( p , t ) ∂ t = g 2 ( x 0 ) Δ x 0 G x 0 ( p , t ) − ∇ V ( x 0 ) · ∇ x 0 G x 0 ( p , t ) − i p U ( x 0 ) G x 0 ( p , t ) , G x 0 ( p , t = 0 ) = 1 , (18) and { ∂ G x 0 ( p , t ) ∂ t = − ( − Δ x 0 ) β / 2 G x 0 ( p , t ) − ∇ V ( x 0 ) · ∇ x 0 G x 0 ( p , t ) − i p U ( x 0 ) G x 0 ( p , t ) . G x 0 ( p , t = 0 ) = 1 . (19) In the following, we provide the forward and the backward Feynman-Kac equations for the microscopic models [ [86] , [90] ]: B ( E α ( t ) ) , B ( S β ( t ) ) , B ( S β ( E α ( t ) ) ) . The forward Feynman-Kac equations: { ∂ G ( x , p , t ) ∂ t = Δ D t 1 − α G ( x , p , t ) − i p U ( x ) G ( x , p , t ) , G ( x , p , t = 0 ) = δ ( x ) , (20) { ∂ G ( x , p , t ) ∂ t = − ( − Δ ) β G ( x , p , t ) − i p U ( x ) G ( x , p , t ) , G ( x , p , t = 0 ) = δ ( x ) , (21) and { ∂ G ( x , p , t ) ∂ t = − ( − Δ ) β D t 1 − α G ( x , p , t ) − i p U ( x ) G ( x , p , t ) , G ( x , p , t = 0 ) = δ ( x ) , (22) which will reduce to the corresponding Fokker-Planck equations when taking p = 0 . The backward Feynman-Kac equations: { ∂ G x 0 ( p , t ) ∂ t = D t 1 − α Δ x 0 G x 0 ( p , t ) − i p U ( x 0 ) G x 0 ( p , t ) , G x 0 ( p , t = 0 ) = 1 , (23) { ∂ G x 0 ( p , t ) ∂ t = − ( − Δ x 0 ) β G x 0 ( p , t ) − i p U ( x 0 ) G x 0 ( p , t ) , G x 0 ( p , t = 0 ) = 1 , (24) and { ∂ G x 0 ( p , t ) ∂ t = D t 1 − α ( − ( − Δ x 0 ) β ) G x 0 ( p , t ) − i p U ( x 0 ) G x 0 ( p , t ) , G x 0 ( p , t = 0 ) = 1 . (25) The forward and backward Feynman-Kac equations for the polymer diffusion with the microscopic model Eq. 5 and Eq. 8 are respectively described as [ [34] , [68] ] { ∂ G ( n , x , p , t ) ∂ t = L n G ( n , x , p , t ) + D ( n ) ∇ x 2 G ( n , x , p , t ) − i p U ( x ) G ( n , x , p , t ) , G ( n , x , p , t = 0 ) = g ( n ) δ ( x ) . (26) and { ∂ G n 0 , x 0 ( p , t ) ∂ t = F n 0 G n 0 , x 0 ( p , t ) + D ( n 0 ) ∇ x 0 2 G n 0 , x 0 ( p , t ) − i p U ( x 0 ) G n 0 , x 0 ( p , t ) , G n 0 , x 0 ( p , t = 0 ) = 1 . (27) 2.2.3. Miscellaneous equations The concept of first passage time, which is the moment when a random variable first attains a specific value, has many applications, e.g., a consideration for investors deciding when to buy or sell stocks during price fluctuations. Escape probability and maximum displacement are closely related concepts to the first passage time and share similar application scenarios. Consider the following stochastic process [ 53 ] { d x ( s ) d s = − ∇ V ( x ( s ) ) + g ( x ( s ) ) ξ 2 ( s ) , d t ( s ) d s = θ ( s ) , (28) where θ ( s ) and ξ 2 ( s ) are independent of each other, t ( s ) is a tempered α -stable Lévy process ( θ ( s ) is the corresponding noise), with the characteristic function given by 〈 e − u t ( s ) 〉 = e − s ( ( u + λ ) α − λ α ) , which results in a α -stable Lévy process when the control factor λ ( ≥ 0 ) approaches 0, and ξ 2 ( s ) is a β -stable Lévy noise. The mean first passage time u ( x ) of the process defined in Eq. 28 from the boundary of Ω starting at the point x ∈ Ω satisfies [ 53 ] { − ∇ V ( x ) · ∇ u ( x ) − ( − Δ ) β / 2 u ( x ) = − α λ α − 1 , x ∈ Ω , u ( x ) = 0 , x ∈ Ω c . (29) Then, we calculate the probability p Γ ( x ) of the jump process defined in Eq. 28 , starting from the point x ∈ Ω , and firstly entering the domain Γ . The p Γ ( x ) is called escape probability solving { − ∇ V ( x ) · ∇ p Γ ( x ) − ( − Δ ) β / 2 p Γ ( x ) = 0 , x ∈ Ω , p Γ ( x ) | x ∈ Γ = 1 , p Γ ( x ) | x ∈ Ω c ∖ Γ = 0 . (30) The microscopic model Eq. 3 can be rephrased as { d v ( t ) d t = − ∇ V ( x ( t ) ) m − γ v ( t ) + ξ 1 ( t ) m , d x ( t ) d t = v ( t ) . (31) The joint distribution f ( x , v , t ) of x ( t ) and v ( t ) solve the equation [ 91 ] ∂ f ( x , v , t ) ∂ t = [ − v ∂ ∂ x + ∂ ∂ v ( γ v + ∇ V ( x ) m ) + 1 m 2 ∂ 2 ∂ v 2 ] f ( x , v , t ) , (32) which is called Klein-Kramers equation. While for the time changed microscopic model v ( E α ( t ) ) and d x ( t ) = d v ( E α ( t ) ) , where v ( t ) is defined as the first equation of Eq. 31 , its Klein-Kramers equation is ∂ f ( x , v , t ) ∂ t + v ∂ ∂ x f ( x , v , t ) = 0 D t 1 − α [ ∂ ∂ v ( γ v + ∇ V ( x ) m ) + 1 m 2 ∂ 2 ∂ v 2 ] f ( x , v , t ) . (33) 3. Analysis and algorithm In this section, we shall analyze the well-posedness of above derived equations in a general form and the regularity of their solutions if possible. Then we propose algorithms, classical numerical algorithms as well as deep learning algorithms, to solve these equations. The analysis of proposed algorithms is performed as well. 3.1. Macroscopic equation for the time-changed strong Markov process All the above discussions start from Brownian motion. Here we discuss a more general form, starting from the strong Markov process X ( t ) . Define the time-changed process Y ( t ) = X ( E ( t ) ) , where E ( t ) is also extended to the inverse of more general subordinator with the characteristic function e − t ϕ ( z ) . It can be noted that when ϕ ( z ) = z α , the subordinator is exactly the above mentioned α -stable Lévy process. The stochastic representation u ( x , t ) = E x [ e − ∫ 0 t κ ( Y ( s ) ) d s f ( Y ( t ) ) ] (34) is the unique mild solution to the equation [ 41 ] { ∂ t w , κ ( x ) u ( x , t ) = L u ( x , t ) − κ ( x ) I t w , κ ( x ) u ( x , t ) , u ( x , 0 ) = f ( x ) , (35) where L is the infinitesimal generator of Y ( t ) . Moreover, the mapping t → u ( · , t ) is continuous with ∥ u ( · , t ) ∥ ≤ C 1 e C 2 t , C 1 , C 2 > 0 , and u ^ ( x , z ) = ∫ 0 ∞ e − z t u ( x , t ) d t , i.e., Laplace transform of u ( x , t ) , exists and satisfies L u ^ ( x , z ) = ϕ ( z + U ( x ) ) u ^ ( x , z ) − f ( x ) ϕ ( z + U ( x ) ) z + U ( x ) . (36) When f ≡ 1 , u x 0 ( p , t ) = E x 0 [ e − i p ∫ 0 t U ( Y ( s ) ) d s ] solves the backward Feynman-Kac equation { ∂ t w , κ ( x 0 ) u x 0 ( p , t ) = L u x 0 ( p , t ) − i p U ( x 0 ) I t w , κ ( x 0 ) u x 0 ( p , t ) , u x 0 ( p , 0 ) = 1 . (37) Applying above results, one can calculate the PDFs of statistical observables (e.g., first passage time, occupation time and path integrals) relating to Y ( t ) . The key here is to get u ( x , t ) from Eq. 35 and Eq. 36 . We take first passage time of non-Markov anomalous subdiffusion as an example and refer interested readers to [ 41 ] for other cases. Precisely, for b ∈ R , define τ b = inf { t > 0 : Y ( t ) > b } , as the considered first passage time. Note that for ρ > 0 , we have P x ( τ b > t ) = P x ( sup 0 ≤ s ≤ t Y ( s ) < b ) = lim ρ → ∞ u ρ ( x , t ) , where u ρ ( x , t ) = E x [ e − ρ ∫ 0 t 1 ( b , ∞ ) ( Y s ) d s ] , is a type of Eq. 34 with κ ( x ) = ρ 1 ( b , ∞ ) ( x ) and f ( x ) ≡ 1 . By subtle techniques and Eq. 36 (see [ 41 ] for details), it is found that u ^ ρ ( x , z ) = 1 z ( 1 − ρ ϕ ( z + ρ ) ( z + ρ ) ( ϕ ( z ) + ϕ ( z + ρ ) ) e ( x − b ) ϕ ( z ) ) . By letting ρ → ∞ and taking inverse Laplace transform, one shall get the probability distribution function of τ b . For many problems, one cannot solve Eq. 35 and Eq. 36 analytically. Hence, one may apply typical numerical methods or fast numerical methods; see [ [92] , [93] , [94] , [95] ] and references therein for constructions of numerical methods and corresponding error estimates. 3.2. Space fractional diffusion equation Consider the two-dimensional space fractional diffusion equation with variable coefficients [ [54] , [55] ] ∂ u ( x , t ) ∂ t = ( d 1 ( x ) x L D x α + d 2 ( x ) x D x R α ) u ( x , t ) + ( e 1 ( x ) y L D y β + e 2 ( x ) y D y R β ) u ( x , t ) + f ( x , t ) , (38) where x : = ( x , y ) , 0 < t ≤ T , α , β ∈ ( 1 , 2 ) are fractional orders, variable coefficients d 1 , d 2 , e 1 , e 2 ≥ 0 , and f are given functions. Here, zero boundary condition and initial condition u ( x , 0 ) = u 0 ( x ) , x ∈ Ω : = ( x L , x R ) × ( y L , y R ) , are imposed. Recall that x L D x α and x D x R α are left and right Riemann-Liouville fractional derivatives [ 96 ] defined by x L D x α u ( x , t ) : = 1 Γ ( 2 − α ) ∂ 2 ∂ x 2 ∫ x L x ( x − ξ ) 1 − α u ( ξ , y , t ) d ξ and x D x R α u ( x , t ) : = 1 Γ ( 2 − α ) ∂ 2 ∂ x 2 ∫ x x R ( ξ − x ) 1 − α u ( ξ , y , t ) d ξ , respectively. Note that y L D y β and y D y R β can be similarly defined. To obtain high order numerical methods for above problems, one of the most important things is to construct high order approximations to left and right Riemann-Liouville fractional derivatives. We have proposed weighted shifted Grünwald-Letnikov discretization (WSGD) [ 55 ] and weighted shifted Lubich discretization (WSLD) [ 54 ] that achieve second order accuracy and fourth order accuracy respectively. Both methods obtain high order accuracy by assembling low order difference operators with appropriate weights and shifts. Since Grünwald-Letnikov difference operator can be regarded as a special case of Lubich difference operator (i.e., fractional linear multistep method [ 92 ]), we may present them in a unified way. The fractional linear multistep method, which is proposed by Lubich in 1986 [ 92 ], can handle Riemann-Liouville type fractional derivatives well with uniform mesh. This method may be characterized by its generating function δ α ( ζ ) = ( ∑ i = 1 L 1 i ( 1 − ζ ) i ) α , where L ≤ 6 and α > 0 . Note that • α = 1 , it reduces to classical L + 1 point backward difference formula; • L = 1 , it is the same as Grünwald-Letnikov difference operator. It is shown in [ 92 ] that the fractional linear multistep method typically gives L -th order accuracy. However, for time dependent problems, a direct application of this method to a space fractional operator with α ∈ ( 1 , 2 ) often results in an unstable numerical scheme. This can be overcomed by shifting the underlying fractional multistep method. Unfortunately, it reduces to first order accuracy for all possible L after shifting. Inspired by linear multistep method, we find that a linear combination of different shifted first order difference operators leads to high order accuracy. More specifically, the shifted Lubich’s difference operator is defined as A h , p α , L u ( x ) : = h − α ∑ k = 0 ∞ q k α , L u ( x − ( k − p ) h ) , where q k α , L is the coefficient of Taylor expansion ( | ζ | ≤ 1 ) δ α ( ζ ) = ∑ k = 0 ∞ q k α , L ζ k . It is only of first order accuracy if p ≠ 0 as least for L ≤ 2 . Assembling the operators with different shifts, high order approximations can be constructed. WSGD method [ 55 ] ( L = 1 ): A linear combination of two different shifts p , q gives L D h , p , q α , 1 u ( x ) = α − 2 q 2 ( p − q ) A h , p α , 1 u ( x ) + 2 p − α 2 ( p − q ) A h , q α , 1 u ( x ) . Under appropriate assumptions, it is proved (mainly using Fourier transform) that this approximation achieves second order accuracy for left RL fractional derivatives. Similar results can be derived for right RL fractional derivatives; see Remark 2.5 in [ 55 ]. Combining this with Crank-Nicolson method gives a high order numerical scheme for space diffusion equation. The numerical scheme is proved to be unconditionally stable when ( p , q ) = ( 1 , 0 ) and ( p , q ) = ( 1 , − 1 ) as the eigenvalues of the matrices corresponding to the discretized operators have negative real parts. Moreover, ∥ u n − U n ∥ ≤ C ( τ 2 + h x 2 + h y 2 ) , where ∥ · ∥ denotes discrete L 2 norm, and τ , h x , h y are stepsizes. As demonstrated in [ 55 ], a linear combination of more than two shifts can achieve at least third order accuracy. However, the resulting numerical scheme may not be unconditionally stable for some α which greatly limits its usage. WSLD method [ 54 ] ( L = 2 ): A linear combination of four different shifts. It can also be obtained by the following way. • second order: a linear combination of two first order approximations ( p ≠ q ) L D h , p , q α , 2 u ( x ) = q p − q A h , p α , 2 u ( x ) + p p − q A h , q α , 2 u ( x ) . • third order: a linear combination of two second order approximations ( r s ≠ p q ) L D h , p , q , r , s α , 2 u ( x ) = w 1 L D h , p , q α , 2 u ( x ) + w 2 L D h , r , s α , 2 u ( x ) , where w 1 = 3 r s + 2 α 3 ( r s − p q ) and w 2 = 3 p q + 2 α 3 ( p q − r s ) . • fourth order: a linear combination of two third order approximations L D h , p , q , r , s , p ¯ , q ¯ , r ¯ , s ¯ α , 2 u ( x ) = w 3 L D h , p , q , r , s α , 2 u ( x ) + w 4 L D h , p ¯ , q ¯ , r ¯ , s ¯ α , 2 u ( x ) , where w 3 and w 4 are constants (see Theorem 2.4 of [ 54 ]). Similar to [ 55 ], a Crank-Nicolson method is used to derive fully discrete numerical scheme. Applying similar techniques as in WSGD method, for some reasonable shifts, the proposed numerical scheme is shown to be unconditionally stable. Moreover, ∥ u n − U n ∥ ≤ C ( τ 2 + h x 4 + h y 4 ) . Following the above constructions, there are two possible extensions to get higher order accuracy: i) starting from fractional multistep method when L ≥ 3 ; ii) using advanced correction techniques to recover high order of shifted operator for L ≥ 2 and then assembling them appropriately. 3.3. Tempered fractional Laplacian From [ 49 ], it is found that the PDF of above mentioned Lévy flight and tempered Lévy flight formally satisfies (see (20)-(22) and (31)-(33) in [ 49 ]) ∂ u ( x , t ) ∂ t = − ( − Δ + λ ) β / 2 u ( x , t ) , x ∈ Ω ⊂ R d , (39) where − ( − Δ + λ ) β / 2 u = c d , β , λ lim ϵ → 0 ∫ R d ∖ B ϵ ( x ) u ( x , t ) − u ( y , t ) e λ | x − y | | x − y | d + β d y is tempered fractional Laplacian of u , with c d , β , λ = Γ ( d / 2 ) 2 π d / 2 | Γ ( − β ) | for λ > 0 and c d , β , 0 = β Γ ( ( d + β ) / 2 ) 2 1 − β π d / 2 Γ ( 1 − β / 2 ) for λ = 0 . There are more discussions relating to tempered fractional Laplacian which is beyond the scope of this work; see [ 97 ] and references therein for details. Note that Eq. 39 is a time dependent problem. Hence, one needs to specify initial and boundary conditions to make it well defined. The initial condition can be easily specified as the value of u ( x , 0 ) in domain Ω . However, specifying the boundary condition is not trivial as the above mentioned Lévy process and tempered Lévy process have discontinuous paths. As a consequence, the majority of trajectories of these stochastic processes cannot hit the boundary ∂ Ω . Therefore, one must account for the information of u ( x , t ) on R d ∖ Ω . In [ 49 ], generalized Dirichlet type and Neumann type boundary conditions are discussed. Here, we focus on the former one. As illustrated in [ 49 ], the appropriate generalized Dirichlet type boundary conditions (refer to (40)-(42) and (49)-(51) in [ 49 ]) are u ( x , t ) | R d ∖ Ω = g ( x , t ) , where for some constants C , M > 0 and small ϵ > 0 , when | x | ≥ M , g ( x , t ) should satisfy | g ( x , t ) | | x | β − ϵ < C , if λ = 0 and | g ( x , t ) | e ( λ − ϵ ) | x | < C , if λ > 0 respectively. Under proper assumptions of g ( x , t ) , the well-posedness of above time dependent equation is proved for λ = 0 (see (71)-(75) in [ 49 ]) and it is claimed that the results for the case λ > 0 can be similarly derived. In some situations, for particles undergoing Lévy flights or tempered Lévy flights, the mean first passage time of particles and escape probability of particles are important. They are closely related to the corresponding steady state equation of the time dependent Eq. 39 which is given by { − ( − Δ + λ ) β / 2 u ( x ) = f ( x ) , x ∈ Ω , u ( x ) = g ( x ) , x ∈ R d ∖ Ω . (40) If f ( x ) = − 1 , g ( x ) = 0 , the solution is the mean first passage time of particles. If f ( x ) = 0 and g ( x ) = { 1 , x ∈ H ⊂ R d ∖ Ω , 0 , x ∈ ( R d ∖ Ω ) ∖ H , the solution represents the probability that particles land in H after first escaping Ω . Now, we discuss the well-posedness of Eq. 40 and the regularity of the solutions of Eq. 40 when g ( x ) = 0 . For simplicity, let us first introduce some notations that is adopted in [ [98] , [99] ]. For 0 < s < 1 , define H s ( Ω ) = { v ∈ L 2 ( Ω ) : | v | H s ( Ω ) < ∞ } , where | v | H s ( Ω ) = ( ∫ ∫ Ω × Ω ( v ( x ) − v ( y ) ) 2 | x − y | d + 2 s d x d y ) 1 / 2 , and H 1 + s ( Ω ) = { v ∈ H 1 ( Ω ) : | ∂ x i v | ∈ H s ( Ω ) , 1 ≤ i ≤ d } . H s ( Ω ) and H 1 + s ( Ω ) are equipped with norm ∥ v ∥ H s ( Ω ) = ∥ v ∥ L 2 ( Ω ) + | v | H s ( Ω ) and ∥ v ∥ H 1 + s ( Ω ) = ∥ v ∥ H 1 ( Ω ) + ∑ i = 1 d | ∂ x i v | H s ( Ω ) respectively. Further, let H s ( Ω ) = { v | Ω : v ∈ H s ( R d ) , v | R d ∖ Ω = 0 } , and H − s ( Ω ) be the dual of H s ( Ω ) . For λ = 0 , by standard techniques, it is shown in [ [98] , [99] ] that for f ∈ H − β / 2 ( Ω ) , there exists a unique solution u ∈ H β / 2 ( Ω ) of the weak formulation of Eq. 40 . Moreover, if f ∈ H r ( Ω ) for some r ≥ − β / 2 , we shall get that u ∈ H β / 2 + γ ( Ω ) with γ = min { β / 2 + r , 1 / 2 − ϵ } for arbitrarily small ϵ > 0 and ∥ u ∥ H β / 2 + γ ( Ω ) ≤ C ∥ f ∥ H r ( Ω ) , holds for smooth domain Ω . If Ω is not smooth (e.g., Lipschitz domain), then u ∈ H β / 2 + 1 / 2 − ϵ ( Ω ) and it holds that ∥ u ∥ H β / 2 + 1 / 2 − ϵ ( Ω ) ≤ C ∥ f ∥ * , where ∥ · ∥ * is some norm (see Theorem 3.3 in [ 99 ]). By defining appropriate weighted Sobolev spaces, an improved regularity result can be obtained (see Theorem 3.5 in [ 99 ]). For λ > 0 , we consider d = 1 [ 56 ]. Suppose that f ∈ H − β / 2 ( Ω ) , the weak formulation of Eq. 40 reads: find u ∈ H β / 2 ( Ω ) such that B ( u , v ) = 〈 f , v 〉 , ∀ v ∈ H β / 2 ( Ω ) , where the duality pair is defined by 〈 f , v 〉 : = ∫ Ω f v d x and the bilinear form is given by B ( u , v ) : = c 1 , β , γ 2 ∫ R ∫ R ( u ( x ) − u ( y ) ) ( v ( x ) − v ( y ) ) e λ | x − y | | x − y | 1 + β d x d y . Similar to [ 98 ], using the Cauchy-Schwarz inequality, it is straightforward to show that the above bilinear form is continuous. Due to the existence of the term e λ | x − y | , the method in [ 98 ] cannot be used directly to prove the coercivity of B ( u , v ) . Fortunately, for Ω = ( a , b ) , we are able to prove that B ( u , v ) is coercive by introducing some new ideas (see Propositions 3.2 and 3.3, Theorem 3.4 in [ 56 ]). Hence, by the Lax-Milgram theorem, Eq. 40 admits a unique solution. Based on the above results, a Riesz basis Galerkin numerical method is proposed in [ 56 ]. Using B-splines M 1 ( x ) and M 2 ( x ) [ 100 ], i.e., M 1 ( x ) = { 1 , x ∈ [ 0 , 1 ] , 0 , otherwise , and M 2 ( x ) = ∫ 0 1 M 1 ( x − t ) d t , the approximation space V n r ( r = 1 or r = 2 ) is constructed as below V n r = { ϕ n , j r ( x ) = 2 n / 2 M r ( 2 n x − j ) , j = 0 , 1 , ⋯ , 2 n − r } , for some n satisfying 2 n ≥ 2 r . We aim to find the numerical solution u n ∈ V n r such that B ( u n , v n ) = 〈 f , v n 〉 , ∀ v n ∈ V n r . For u ∈ H μ ( Ω ) ∩ H β / 2 ( Ω ) ( μ ≥ β / 2 ), it is shown that the numerical solution has the following approximation property ∥ u − u n ∥ H β / 2 ( R ) ≤ C 2 − n ( min μ , r − β / 2 ) ∥ u ∥ H β / 2 ( Ω ) . 3.4. Tempered fractional Feynman-Kac equation Analyzing functional distribution of tempered anomalous dynamics is one of the feasible approaches to characterize it. As illustrated in Section 2 , such a functional distribution is typically governed by the tempered fractional Feynman-Kac equation [ [51] , [57] ]. In this part, we aim to introduce an efficient time discretization method for solving the following equation D t ( x ) u ( x , t ) − ( λ α + Δ ) D t 1 − α u ( x , t ) = − u 0 ( x ) ( λ α D t 1 − α − λ ) e i p U ( x ) t (41) with initial condition u ( x , 0 ) = u 0 ( x ) , x ∈ Ω ⊂ R d , and zero boundary condition. Here, p is the parameter representing the characteristic function of the joint probability of ( x ( t ) , A ) with A = ∫ 0 t U ( x ( τ ) ) d τ , D t ( x ) is the substantial derivative given by D t ( x ) : = λ − i p U ( x ) + ∂ ∂ t , and D t 1 − α ( 0 < α < 1 ) is the Riemann-Liouville fractional substantial derivative given by D t 1 − α u = 1 Γ ( α ) D t ( x ) ∫ 0 t e − ( t − s ) ( λ − i p U ( x ) ) ( t − s ) 1 − α u ( x , s ) d s . Let us start with the well-posedness of Eq. 41 . The analysis is based on the integral representation of u , which is derived mainly using Laplace transform and inverse Laplace transform. In fact, taking Laplace transform of Eq. 41 gives ( η ( z ) − Δ ) β ( z ) 1 − α u ^ ( x , z ) = u 0 ( x ) β ( z ) 1 − α η ( z ) z − i p U ( x ) , (42) where u ^ ( x , z ) is the Laplace transform of u ( x , t ) , β ( z ) , and η ( z ) are given by β ( z ) = z + λ − i p U ( x ) and η ( z ) = β ( z ) α − λ α , respectively. After getting an explicit expression of u ^ ( x , z ) and then taking inverse Laplace transform, one obtains u ( x , t ) = 1 2 π i ∫ Γ θ , κ e z t β ( z ) α − 1 ( η ( z ) − Δ ) − 1 ( β ( z ) 1 − α u 0 ( x ) η ( z ) z − i p U ( x ) ) d z , where Γ θ , κ is a contour (see (2.25) in [ 57 ]). Under appropriate assumptions, it is proved that the above integral representation is the mild solution of tempered fractional Feynman-Kac Eq. 41 (see Proposition 3.1(3) in [ 57 ]). One can note that the analysis is one of the main results in [ 57 ] and the integral expression is very important in the analysis of the numerical method in what follows. Now we present the time-stepping scheme that is derived from the discretization in frequency domain, i.e., discretization based on Eq. 42 . To this aim, we first construct an approximation of D t 1 − α u ( x , t n ) in spatial domain. Then the final discretization is motivated by relating a transformation of this approximation to the Laplace transform of D t 1 − α u ( x , t ) . More specifically, we have • Approximation of D t 1 − α u ( x , t n ) . By straightforward calculations, it is observed that D t 1 − α u ( x , t ) = e − t ( λ − i p U ( x ) ) ( ∂ ∂ t ) 1 − α ( e t ( λ − i p U ( x ) ) u ( x , t ) ) , where ( ∂ ∂ t ) 1 − α u is the standard Riemann-Liouville fractional derivative defined by ( ∂ ∂ t ) 1 − α u ( x , t ) = 1 Γ ( α ) ∂ ∂ t ∫ 0 t ( t − s ) α − 1 u ( x , s ) d s . Approximating ( ∂ ∂ t ) 1 − α u by backward Euler convolution quadrature ∂ ¯ τ 1 − α u n = 1 τ 1 − α ∑ j = 1 n b n − j ( 1 − α ) u j , where τ is the time step size and b n − j ( 1 − α ) is the coefficient (see (2.5) in [ 57 ]), one can obtain an approximation of D t 1 − α u ( x , t n ) , that is, D ¯ τ 1 − α u n ( x ) : = e − t n ( λ − i p U ( x ) ) ∂ ¯ τ 1 − α ( e t n ( λ − i p U ( x ) ) u ( x , t ) ) . • Transformation of D ¯ τ 1 − α u n ( x ) . It is found that ∑ n = 1 ∞ D ¯ τ 1 − α u n ( x ) ζ n = ( 1 − e − τ ( λ − i p U ( x ) ) ζ τ ) 1 − α ∑ n = 1 ∞ u n ( x ) ζ n . Noticing that the Laplace transform of D t 1 − α u ( x , t ) is given by β ( z ) 1 − α u ^ ( x , z ) = ∫ 0 ∞ D t 1 − α u ( x , t ) e − t z d t ≈ τ ∑ n = 1 ∞ D ¯ τ 1 − α u n ( x ) e − t n z . Taking ζ = e − τ z and comparing above results motivate us to consider the following approximations β ( z ) ≈ 1 − e − τ ( z + λ − i p U ( x ) ) τ , u ^ ( x , z ) ≈ τ ∑ n = 1 ∞ u n ( x ) e − t n z . We also apply above approximation of β ( z ) in η ( z ) . For 1 / ( z − i p U ( x ) ) , we use τ e − τ ( z − i p U ( x ) ) / ( 1 − e − τ ( z − i p U ( x ) ) ) instead of τ / ( 1 − e − τ ( z − i p U ( x ) ) ) for analysis purpose. Based on the above discussions, the resulting numerical scheme is obtained, that is, ( D ¯ τ α − λ α − Δ ) D ¯ τ 1 − α u n ( x ) = u 0 D ¯ τ 1 − α ( D ¯ τ α − λ α ) e i p U ( x ) t n , n ≥ 1 . Applying Cauchy’s integral formula, we are able to get an integral expression of u n ( x ) . With the help of explicit expressions of u ( x , t ) and u n ( x ) , it is proved that (see Theorem 3.3 in [ 57 ]) ∥ u ( · , t n ) − u n ∥ M ( Ω ) ≤ C ∥ u 0 ∥ M ( Ω ) t n − 1 τ , n ≥ 1 , where ∥ · ∥ M ( Ω ) denotes the dual norm of C ( Ω ¯ ) . 3.5. Equation driven by fractional Gaussion noise Let us first briefly describe the equation considered in this part. Assume that Ω is a bounded domain, let B ( t ) be a standard Brownian motion with B ( 0 ) ∈ Ω and S ( t ) be a β -stable subordinator. Define a stochastic process X ( t ) as follows X ( t ) = { B ( S ( t ) ) , S ( t ) ≤ τ Ω , Θ , S ( t ) ≥ τ Ω , where τ Ω = inf { t > 0 : B ( t ) ∉ Ω } is a stopping time of B ( t ) and Θ is a coffin state. From [ 101 ], it is seen that the infinitesimal generator of X ( t ) is the spectral fractional Laplacian operator ( − Δ ) β ( β ∈ ( 0 , 1 ) ) that is defined by ( − Δ ) β u = ∑ k = 1 ∞ λ k β ( u , ϕ k ) ϕ k , where { λ k , ϕ k } k = 1 ∞ are eigenvalues and eigenfunctions (in L 2 ( Ω ) ) pairs of − Δ with zero boundary condition. Then, the Fokker-Planck equation relating to X ( t ) time changed by the inverse α -stable subordinator is ∂ u ∂ t + ( ∂ ∂ t ) 1 − α ( − Δ ) β u = 0 . If there exists external fractional Gaussian noise and external source term depending on the density of particles, then we get the fractional diffusion equation driven by fractional Gaussian noise [ 58 ] ∂ u ∂ t + ( ∂ ∂ t ) 1 − α ( − Δ ) s u = f ( u ) + W Q H ˙ , x ∈ Ω , t ∈ ( 0 , T ] (43) with zero initial and boundary conditions. Here, f is a nonlinear term satisfying ∥ f ( u ) ∥ L 2 ( Ω ) ≤ C ( 1 + ∥ u ∥ L 2 ( Ω ) ) , ∥ f ( u ) − f ( v ) ∥ L 2 ( Ω ) ≤ C ∥ u − v ∥ L 2 ( Ω ) . W Q H is the fractional Gaussian process defined by W Q H = ∑ k = 1 ∞ Λ k ϕ k W k H , where { W k H } k = 1 ∞ are one dimensional fractional Brownian motions that are mutually independent, H ∈ ( 0 , 1 ) is the Hurst index, and Q is a nonnegative linear self-adjoint operator that has the same eigenfunctions with − Δ . The corresponding eigenvalues of Q are denoted by { Λ k } k = 1 ∞ . Now we are ready to discuss the regularity of the mild solution of Eq. 43 . First, an equivalent expression of the mild solution is given by taking Laplace transform and inverse Laplace transform, that is, u = ∫ 0 t R ( t − r ) f ( u ( r ) ) d r + ∫ 0 t R ( t − r ) d W Q H ( r ) = ∫ 0 t R ( t − r ) f ( u ( r ) ) d r + ∑ k = 1 ∞ ∫ 0 t Λ k E k ( t − r ) ϕ k d W k H ( r ) , where R ( t ) = 1 2 π i ∫ Γ θ , κ e z t z α − 1 ( z α + ( − Δ ) s ) − 1 d z and E k ( t ) = 1 2 π i ∫ Γ θ , κ e z t z α − 1 ( z α + λ k s ) − 1 d z with Γ θ , κ being a contour as above. Then applying estimates of R ( t ) and E k ( t ) (see (2.5) and (2.6) in [ 58 ]) and the regularization of noise (see Lemma 2.6 in [ 58 ]) E [ ( ∫ 0 T g ( T − r ) d W k H ( r ) ) 2 ] ≤ C ∥ ( ∂ ∂ t ) 1 / 2 − H g ∥ L 2 ( [ 0 , T ] ) 2 , one can get the spatial regularity and temporal Hölder regularity of the mild solution u under appropriate assumptions, i.e., E [ ∥ ( − Δ ) σ u ∥ L 2 ( Ω ) 2 ] ≤ C , E [ ∥ u ( t ) − u ( t − τ ) τ γ ∥ L 2 ( Ω ) 2 ] ≤ C , for some σ and γ (see Theorems 2.8 and 2.9 in [ 58 ]). Next, we present the numerical method in [ 58 ]. The main idea is to apply spectral Galerkin method to discretize the fractional Laplacian and backward Euler convolution quadrature to discretize ( ∂ ∂ t ) 1 − α u . The procedure is as follows: 1. Semidiscrete scheme. Using the above eigenfunctions, define a finite dimensional space V N as V N = span { ϕ 1 , ⋯ , ϕ N } ⊂ L 2 ( Ω ) . Our objective is to find u N ( t ) ∈ V N such that { ∂ u N ∂ t + ( ∂ ∂ t ) 1 − α ( − Δ ) N s u N = P N f ( u N ) + P N W Q H ˙ , u N ( 0 ) = 0 , (44) where P N u = ∑ i = 1 N ( u , ϕ i ) ϕ i is a projection of L 2 ( Ω ) onto V N and ( − Δ ) N s : V N → V N is defined by ( ( − Δ ) N s u N , v N ) = ( ( − Δ ) s u N , v N ) , ∀ v N ∈ V N . Obviously, Eq. 44 has a similar form to Eq. 43 . Hence, by similar techniques, one can obtain an explicit expression of u N ( t ) as well as the corresponding estimates. Further, one can get the spatial error E [ ∥ u − u N ∥ L 2 ( Ω ) 2 ] 1 / 2 ≤ C ( N + 1 ) − 2 σ / d . 2. Fully discrete scheme. The further discretization is based on Eq. 44 . In fact, applying standard finite difference to ∂ u ∂ t and W Q H ˙ in Eq. 44 , and backward Euler method to time fractional term ( ∂ ∂ t ) 1 − α in Eq. 44 , one can get the fully discrete scheme of Eq. 43 ∂ τ u N n + ∂ ¯ τ 1 − α ( − Δ ) N s u N n = P N f ( u N n − 1 ) + P N ∂ τ W Q H ( t n ) , where ∂ τ u ( t n ) = u ( t n ) − u ( t n − 1 ) τ and ∂ ¯ τ 1 − α is the same as before. With the help of a proper transformation, one can also obtain an explicit expression of u N n in a form similar to the one of u N ( t ) . In this manner, for sufficiently small ϵ > 0 , one can derive the following temporal error E [ ∥ u N ( t n ) − u N n ∥ L 2 ( Ω ) 2 ] 1 / 2 ≤ C τ H − ρ α / s − ϵ , where 0 < ρ < min { s H / α , s } with 0 < α < 1 . Using the above derived spatial and temporal error, by triangle inequality, the numerical approximation error is given by E [ ∥ u ( t n ) − u N n ∥ L 2 ( Ω ) 2 ] 1 / 2 ≤ C ( ( N + 1 ) − 2 σ / d + τ H − ρ α / s − ϵ ) . 3.6. Semilinear parabolic PDEs with infinite dimensional coupling In [ 68 ], Fokker-Planck equations and Feynman-Kac equations that describe some statistical observables of polymer dynamics models are derived. These equations can be reformulated as semilinear parabolic PDEs with infinite dimensional coupling { ∂ u ( n , x , t ) ∂ t + T n u ( n , x , t ) + T x u ( n , x , t ) + f = 0 , u ( n , x , T ) = g ( n , x ) . (45) Here u : N × R d × [ 0 , ∞ ) → R is unknown, T n and T x are operators defined by T n u ( n ) : = { α ( 0 ) ( u ( 1 ) − u ( 0 ) ) , n = 0 , α ( n ) ( u ( n + 1 ) − u ( n ) ) + β ( n ) ( u ( n − 1 ) − u ( n ) ) , n ≥ 1 , where α ( n ) , β ( n ) are given functions and T x : = 1 2 Tr ( ( σ σ T ) ( n , x , t ) H e s s x ) + μ ( n , x , t ) · ∇ x with μ ∈ R d , σ ∈ R d × d being vector-valued and matrix-valued functions, respectively, g and f = f ( t , n , x , u ( n , x , t ) , σ T ∇ x u ( n , x , t ) ) being scalar-valued functions. Because of the operator T n , Eq. 45 becomes an infinite dimensional coupled lattice system. It is extremely difficult to solve it by traditional numerical methods if it is not impossible. Hence, we turn to deep neural network method. To be specific, we are going to extend standard deep BSDE method [ [64] , [65] ] that works well for high dimensional nonlinear PDEs to infinite dimensional systems Eq. 45 . If there is no operator T n , then standard deep BSDE method can be applied directly to Eq. 45 . Roughly speaking, one needs three steps to construct a standard deep BSDE method [ 65 ], i.e., • Step 1: construct a stochastic process X ( t ) , d X ( t ) = μ d t + σ d B ( t ) , X ( 0 ) = x , the infinitesimal generator of which is exactly T x . • Step 2: derive a BSDE by applying Itô’s formula to u ( X ( t ) , t ) first and then replacing u ( X ( t ) , t ) and σ ∇ x u ( X ( t ) , t ) by new notations Y ( t ) and Z ( t ) . In this manner, one can get d Y ( t ) = − f ( t , X ( t ) , Y ( t ) , Z ( t ) ) d t + Z ( t ) T d B ( t ) , Y ( T ) = g ( X ( T ) ) with ( Y ( t ) , Z ( t ) ) = ( u ( X ( t ) , t ) , σ ∇ x u ( X ( t ) , t ) ) being the unique solution. This together with Step 1 gives so-called forward-backward stochastic differential equation (FBSDE) { X ( t ) = X ( 0 ) + ∫ 0 t μ d s + ∫ 0 t σ d B ( s ) , Y ( t ) = Y ( T ) + ∫ t T f d s − ∫ t T Z ( s ) T d B ( s ) X ( 0 ) = x , Y ( T ) = g ( X ( T ) ) , Hence, Eq. 45 without operator T n can be formulated as a constrained optimization problem below inf Y 0 , { Z t } 0 ≤ t ≤ T E [ | Y T − g ( X T ) | 2 ] , such that { X ( t ) = X ( 0 ) + ∫ 0 t μ d s + ∫ 0 t σ d B ( s ) , Y ( t ) = Y ( 0 ) − ∫ 0 t f d s + ∫ 0 t Z ( s ) T d B ( s ) . • Step 3: solve above optimization problem by discretizing constraints and approximating Y ( 0 ) and Z ( t ) via independent neural networks. Since Y ( 0 ) = u ( X ( 0 ) , t = 0 ) = u ( x , 0 ) , one can get an approximation of u ( x , 0 ) . The deep learning method in [ 68 ] follows a similar procedure as above and handles issues incurred by operator T n . The first issue is to find an appropriate stochastic process whose infinitesimal generator is T n . Unfortunately, it is not easy to find such a process. This may be overcomed partially by its microscopic description, i.e., the birth-death process N ( t ) that satisfies Eq. 8 . Since N ( t ) is essentially a jump process, we are not able to eliminate T n totally as T x . In fact, applying Itô’s formula to u ( N ( t ) , X ( t ) , t ) , one can get the following BSDE (see Lemma 2.1 and Theorem 2.1 in [ 68 ]) d u ( N ( t ) , X ( t ) , t ) = − f d t + [ ∇ x u ] T σ d B ( t ) + ∫ Z ∖ { 0 } δ u ( t , n ; N ( t − ) ) J ˜ ( d t , d n ; N ( t − ) ) (46) with terminal condition u ( N ( T ) , X ( T ) , T ) = g ( N ( T ) , X ( T ) ) . Here δ u ( t , n ; N ( t − ) ) = u ( N ( t − ) + n , X ( t ) , t ) − u ( N ( t − ) , X ( t ) , t ) and J ˜ ( d t , d n ; N ( t − ) ) is a compensated counting random measure (see (2.4)-(2.6) in [ 68 ]). Obviously, one can approximate u ( N ( 0 ) , X ( 0 ) , 0 ) and [ ∇ x u ] T σ in Eq. 46 by independent neural networks as in Step 3. The remaining term in Eq. 46 is the second issue that we need to deal with. It can be noted that for the above birth-death process ∫ Z ∖ { 0 } δ u ( t , n ; N ( t − ) ) J ˜ ( d t , d n ; N ( t − ) ) = [ u ( N ( t ) , X ( t ) , t ) − u ( N ( t − ) , X ( t ) , t ) ] d t − α ( N ( t − ) ) [ u ( N ( t − ) + 1 , X ( t ) , t ) − u ( N ( t − ) , X ( t ) , t ) ] d t − β ( N ( t − ) ) [ u ( N ( t − ) − 1 , X t , t ) − u ( N ( t − ) , X ( t ) , t ) ] d t , where u ( N ( t − ) ± 1 , X ( t ) , t ) are unknown. One may also approximate them using neural networks as in Step 3. For simplicity, we devise a vector-valued neural network to approximate δ u ( t , ± 1 ; N ( t − ) ) directly. Till now, we have addressed two main issues encountered when applying deep BSDE method to Eq. 45 . The full deep BSDE method for Eq. 45 can be derived easily following steps 1–3. We omit details here and refer interested readers to our work [ 68 ]. 4. Applications in chemistry and biology 4.1. Modeling telomere shortening process Aging is a complex biological process influenced by genes, environment, and lifestyle; and investigating its molecular and cellular changes can uncover potential mechanisms and intervention strategies. Telomeres are specific DNA sequences at the ends of linear chromosomes, and their shortening is associated with cellular aging, death, and cancer. However, some cells combat this process by expressing telomerase to repair and lengthen telomeres [ [102] , [103] ]. Telomere shortening (TS) is primarily caused by incomplete replication of chromosomes, the action of exonucleases, and damage induced by oxidative stress; this application of above discussions will simulate the dynamic behavior of telomere length, derive macroscopic equations, and calculate the distribution of relevant statistical measures [ 104 ]. In the medical field, the measurement of telomere length serves as a crucial diagnostic tool [ 105 ]. To comprehend the complex processes of TS at the microscopic level, researchers have examined the underlying mechanisms of telomere length dynamics from a stochastic perspective. According to [ 7 ], incomplete replication of chromosome ends leads to TS, with the shortened length L 1 following a normal distribution. The probability of TS due to exonuclease activity is measured/assumed to be 0.8, and the shortened length L 2 follows a Poisson distribution. During cellular replication, oxidative stress causes DNA damage to telomeres, leading to a measured/assumed probability of TS of 0.1, and the shortened length L 3 also follows a normal distribution. Then, there holds L = 1 × 1 4 × L 1 + 0.8 × 1 4 × L 2 + 0.1 × 1 2 × L 3 × N , (47) where L represents the shortened length of the telomere and N is the number of bases damaged in the DNA strand. It is assumed that the waiting time adheres to a tempered power-law distribution, and that the TS jump length L changes independently at each step. Considering the low probability of TS due to oxidative stress damage, only the effects of incomplete replication at chromosome ends and exonuclease activity are taken into account. Therefore, L = L 1 + L 2 , and it is assumed that L 1 and L 2 are independent. Then, one can get the PDF of L as φ ( L ) = φ 1 ( L 1 ) * φ 2 ( L 2 ) ∼ ∑ L 2 = 0 ∞ [ 1 2 π σ exp { − ( L − μ − L 2 ) 2 2 σ 2 } × exp { − η } η L 2 L 2 ! ] , (48) where “ * ” represents the convolution operation, and η is the intensity of the Poisson distribution, μ and σ are the mean and variance of the normal distribution, respectively. Considering that the functional of length L is the conditional probability density G of L ( 0 ) = L 0 , one can derive the backward Feynman-Kac equation [ 104 ] ∂ G L 0 ( p , t ) ∂ t = D t 1 − α , λ [ ( η + μ ) ( 1 + B α λ α ) B α ∂ ∂ L 0 + σ 2 2 B α ∂ 2 ∂ L 0 2 + λ α ] G L 0 ( p , t ) − [ λ + p U ( L 0 ) ] G L 0 ( p , t ) + ( λ − λ α D t 1 − α , λ ) e − p U ( L 0 ) t . (49) Since telomere length is not infinite at the onset of a cell’s life but starts with an initial length of l 0 . When telomeres shorten to a certain degree, the stability of the genome within the cell is compromised, ultimately leading to cellular aging, death, or cancer. Specifically, when the length of the shortest telomere in the cell reaches the critical threshold l c , the cell’s capacity for division becomes restricted and it begins to senesce. Therefore, the upper bound of the shortened length is l 0 − l c . The occupation time is the total time for a telomere to shorten the length between [ 0 , l 0 − l c ] in the observation time [ 0 , t ] , which can be defined as [ 104 ] T + = ∫ 0 t U ( L ( τ ) ) τ , (50) where U ( L ) = { 1 , L ∈ [ 0 , l 0 − l c ] , 0 , L ∉ [ 0 , l 0 − l c ] . (51) Exploring the total time of TS deepens our insight into cell aging mechanisms. TS is closely tied to the development of diseases like cancer, cardiovascular issues, and neurological conditions, offering key clues for prevention and treatment strategies. This application aids in anti-aging research and understanding disease progression. In [ 104 ], Figure 4 shows J ( t ) peaking before declining to zero, indicating the time most cell telomeres reach l c . Due to the monotonically decreasing distribution of TS jump lengths, the occupancy time shares the same shape as the distribution of the first passage time. 4.2. Time-changed tempered fractional Langevin-Brownian motion In certain real-world datasets, such as those in biology [ [106] , [107] ], financial time series [ 108 ], ecology [ 109 ], and physics [ 110 ], a time-changed stochastic process is required. This process involves substituting the deterministic time variable with a positive, non-decreasing random process, which results in a blend of two independent random processes. One of these processes is referred to as the external process (or the original process), while the other is known as the internal process (or a subordinator). Tempered fractional Langevin equation is driven by tempered fractional Gaussian noise γ ( t ) [ 44 ]. It is also a Gaussian process and can be written as { d x ( s ) d s = v ( s ) , d v ( s ) d s = − ∫ 0 s K ( s − τ ) v ( τ ) d τ + ρ γ ( s ) , d t ( s ) d s = η ( s ) , (52) where ρ = 2 k B T , the kernel K ( t ) = 2 〈 γ ( 0 ) γ ( t ) 〉 = h − 2 ( C t + h 2 | t + h | 2 H + C t − h 2 | t − h | 2 H − 2 C t 2 | t | 2 H ) for a sufficient small h , C t 2 = 2 Γ ( 2 H ) ( 2 λ | t | ) 2 H − 2 Γ ( H + 1 2 ) K H ( λ | t | ) π ( 2 λ | t | ) H , and K H ( t ) is the modified Bessel function of second kind. We assume the initial velocity satisfies the condition v 0 2 = k B T . The PDF of the subordinated process X ( t ) : = x ( s ( t ) ) can be written as p ( x , t ) = ∫ 0 ∞ p 0 ( x , s ) f ( s , t ) d s , (53) where p 0 ( x , s ) is the PDF of the original process x ( s ) and f ( s , t ) is the PDF of the inverse β -stable subordinator s ( t ) . The moments of subordinated process X ( t ) could be obtained by the relation L t → u 〈 X n ( t ) 〉 = u β − 1 L s → u β 〈 x n ( s ) 〉 (54) in Laplace space. According to Eq. 54 , with the time evolution the first and second moments of the subordinated process X ( t ) : = x ( s ( t ) ) behave as 〈 X ( t ) 〉 : k B T β Γ ( β ) t β (55) and 〈 X 2 ( t ) 〉 : k B T β Γ ( 2 β ) t 2 β → F t ( 2 − 2 H ) β → k B T A β Γ ( 2 β ) t 2 β , (56) where E = k B T / [ 2 D H Γ 2 ( H + 1 / 2 ) Γ ( 2 H + 1 ) Γ ( ( 1 − 2 H ) β + 1 ) ] and F = k B T / [ D H Γ 2 ( H + 1 / 2 ) Γ ( 2 H + 1 ) Γ ( ( 2 − 2 H ) β + 1 ) ] . The MSD of the time-changed tempered fractional Langevin equation evolves over time as [ 46 ] 〈 ( Δ X ( t ) ) 2 〉 : ( k B T β Γ ( 2 β ) − k B T β 2 Γ 2 ( β ) ) t 2 β → F t ( 2 − 2 H ) β − E 2 t 2 ( 1 − 2 H ) β → ( k B T A β Γ ( 2 β ) − A 2 ( β Γ ( β ) ) 2 ) t 2 β . (57) This implies that the time-changed tempered fractional Langevin process exhibits different characteristics at different time scales. 5. Future prospects 5.1. Biological macromolecules dynamics The four most essential macromolecules in living organisms are nucleic acids, proteins, polysaccharides, and lipids. Materials like plastics, rubber, and fibers, which are types of polymer materials, have dramatically changed our daily life. These substances are created through the processes of polymerization and depolymerization from similar monomers. When studying the kinetic behavior of these materials, it’s crucial to take into account not only the polymer’s inherent movement characteristics but also the effects of chemical interactions during polymerization and depolymerization, as well as the influence of the surrounding environment. To delve deeper into the kinetic behavior of cell division and intact polymer proteins within living organisms, the following research initiatives are planned for future exploration. 5.1.1. Kinetic modeling of microtubules The dynamics and control of microtubules are vital for the proper functioning and division of all eukaryotic cells [ [111] , [112] ]. As depicted in [ 111 ], microtubules extend and attach to the replicated DNA, forming a spindle and generating the pulling force that initiates cell division. Therefore, investigating the growth and regulation at the ends of microtubules can enhance our understanding of the mitotic behavior in eukaryotic cells. 5.1.2. Protein synthesis, transport, and movement Protein serves as a fundamental polymer in the human body, and its synthesis is a direct outcome of gene expression, which encompasses a series of processes such as transcription and translation. During the gene transcription process, RNA polymerase II moves along the DNA sequence, synthesizing mRNA. This movement can exhibit three distinct phases: transcriptional elongation, backtracking, and regressive recovery [ 113 ]. The translation process takes place in the ribosomes, hence requiring the transportation of the transcribed mRNA from the nucleus to the intended location following the completion of transcription [ 4 ]. The diffusion of synthesized proteins, like neurotransmitter receptors, on the surface of cell membranes can be influenced by crowded environments (as depicted in [ 5 ]). Therefore, it is both interesting and important to consider the modeling of the periodic behavior of proteins. 5.2. Multi-fluid modeling 5.2.1. Modeling sediment transport by wind The collective process of sand and dust being emitted, transported, and deposited by the wind is known as aeolian processes, named after the Greek god Aeolus, who was the keeper of the winds [ [114] , [115] ]. Aeolian processes occur in areas where there is an ample supply of granular material and winds of sufficient force to move them through the atmosphere. On Earth, this phenomenon is most pronounced in deserts, on beaches, and in other areas with sparse vegetation, such as dried-up lake beds. The blowing of sand and dust in these regions plays a pivotal role in shaping the landscape through the formation of sand dunes and ripples, the erosion of rocks, and the creation and transport of soil particles. Furthermore, airborne dust particles can be carried for thousands of kilometers from their original source, impacting weather and climate, ecosystem productivity, the hydrological cycle, and various other components of the Earth’s system. Consequently, the study of the kinetic behavior of wind-blown sand is a significant area of research. 5.2.2. Modeling fluid and solid interaction The study of microclimates within plant canopies has long been a source of inspiration for scientists engaged in diverse research fields, including agronomy, ecology, and silviculture. It was nearly a century ago that the first measurements of wind speed within a forest stand were published in [ 116 ]. The behavior of wind in the canopy is an important component of the canopy microclimate, which largely determines the rate of exchange of heat, water vapor, and other associated gases and particles with the atmosphere. Consequently, the second topic of study is the interactions in the canopy and the wind field, which can be of great assistance in wind and sand control, seed dispersal [ 117 ], and in understanding inversions in agriculture. 5.2.3. Modeling wind and fluid interaction, and the aroused enhanced diffusion Ocean-atmosphere interactions exert a significant influence on the marine environment. For instance, hurricanes can impact upper ocean temperatures [ [118] , [119] ], while interactions between ocean currents and winds affect surface carbon concentrations and air-sea carbon exchange in the Southern Ocean [ 120 ]. Global warming can be interrupted by the Pacific circulation [ 121 ], and the interactions between the ocean and wind can directly impact the dispersion of marine pollutants [ 122 ] and more. Therefore, general coupled ocean-atmosphere and pollutant dispersion modeling is an important and intriguing research topic. 5.3. General form of chemotaxis model The myxobacteria are ubiquitous soil bacteria that aggregate under conditions of starvation and construct fruiting bodies as a means of survival. The mechanisms underlying their social gliding, aggregation, and fruiting body formation have remained poorly understood until recently. In [ 123 ], a stochastic cellular automaton model is presented with the objective of describing and providing an understanding of the mechanisms by which the bacteria manage to build higher-organized structures. This model is affected by three factors which are, slime, diffusing chemoattractant, and inertia of motion. Therefore, an interesting topic is to derive the equations satisfied by the statistical observables of this chemotaxis model. By studying the equations, one can well understand the chemotaxis phenomena. Declaration of competing interest The authors declare that they have no conflicts of interest in this work. Acknowledgments This work was supported by the National Natural Science Foundation of China (12225107 and 12071195), the Major Science and Technology Projects in Gansu Province-Leading Talents in Science and Technology (23ZDKA0005), the Innovative Groups of Basic Research in Gansu Province (22JR5RA391), and Lanzhou Talent Work Special Fund. Biographies Heng Wang, PhD Student of School of Mathematics and Statistics, Lanzhou University. His research interests include machine learning, AI for science, data and mechanism integration algorithm. Weihua Deng ( BRID:09663.00.71823 ), Professor of School of Mathematics and Statistics and National Key Laboratory of Natural Product Chemistry, Lanzhou University. His research interests include multi-scale modeling, scientific computation, and deep learning. Footnotes Peer review under the responsibility of Editorial Board of Fundamental Research. References 1. Wang D., Bushnell D.A., Huang X., et al. Structural basis of transcription: Backtracked RNA polymerase II at 3.4 angstrom resolution. Science. 2009;324(5931):1203–1206. doi: 10.1126/science.1168729. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Cheung A., Cramer P. Structural basis of RNA polymerase II backtracking, arrest and reactivation. Nature. 2011;471:249–253. doi: 10.1038/nature09785. [ DOI ] [ PubMed ] [ Google Scholar ] 3. Gonzalez M.N., Blears D., Svejstrup J.Q. Causes and consequences of RNA polymerase II stalling during transcript elongation. Nat. Rev. Mol. Cell Biol. 2021;22:3–21. doi: 10.1038/s41580-020-00308-8. [ DOI ] [ PubMed ] [ Google Scholar ] 4. Jansen R. mRNA localization: Message on the move. Nat. Rev. Mol. Cell Biol. 2001;2:247–256. doi: 10.1038/35067016. [ DOI ] [ PubMed ] [ Google Scholar ] 5. Mosqueira A., Camino P.A., Barrantes F.J. Antibody-induced crosslinking and cholesterol-sensitive, anomalous diffusion of nicotinic acetylcholine receptors. J. Neurochem. 2020;152(6):663–674. doi: 10.1111/jnc.14905. [ DOI ] [ PubMed ] [ Google Scholar ] 6. Ye Z., Zhang C., Yuan J., et al. Ligand-receptor interaction triggers hopping and sliding motions on living cell membranes. J. Am. Chem. Soc. 2023;145(46):25177–25185. doi: 10.1021/jacs.3c06925. [ DOI ] [ PubMed ] [ Google Scholar ] 7. Zhang Q.H., Tian X.J., Liu F., et al. A switch-like dynamic mechanism for the initiation of replicative senescence. FEBS Lett. 2014;588(23):4369–4374. doi: 10.1016/j.febslet.2014.09.043. [ DOI ] [ PubMed ] [ Google Scholar ] 8. Ishikawa N., Nakamura K.-I., Izumiyama-Shimomura N., et al. Changes of telomere status with aging: An update. Geriatr. Gerontol. Int. 2016;16(S1):30–42. doi: 10.1111/ggi.12772. [ DOI ] [ PubMed ] [ Google Scholar ] 9. Baldovin F., Orlandini E., F S. Polymerization induces non-Gaussian diffusion. Front. Phys. 2019;7(124) [ Google Scholar ] 10. Banks D.S., Fradin C. Anomalous diffusion of proteins due to molecular crowding. Biophys. J. 2005;89(5):2960–2971. doi: 10.1529/biophysj.104.051078. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Hu Q., Ewing R.P., Dultz S. Low pore connectivity in natural rock. J. Contam. Hydrol. 2012;133:76–83. doi: 10.1016/j.jconhyd.2012.03.006. [ DOI ] [ PubMed ] [ Google Scholar ] 12. Bartumeus F., da Luz M.G.E., Viswanathan G.M., et al. Animal search strategies: A quantitative random-walk analysis. Ecology. 2005;86(11):3078–3087. [ Google Scholar ] 13. Brockmann D., Hufnagel L., Geisel T. The scaling laws of human travel. Nature. 2006;439:462–465. doi: 10.1038/nature04292. [ DOI ] [ PubMed ] [ Google Scholar ] 14. Montroll E.W., Weiss G.H. Random walks on lattices. II. J. Math. Phys. 1965;6(2):167–181. [ Google Scholar ] 15. Berkowitz B., Cortis A., Dentz M., et al. Modeling non-Fickian transport in geological formations as a continuous time random walk. Rev. Geophys. 2006;44(2) [ Google Scholar ] 16. Zhang J., Chen A., Qiu H., et al. Exact results for gene-expression models with general waiting-time distributions. Phys. Rev. E. 2024;109:024119. doi: 10.1103/PhysRevE.109.024119. [ DOI ] [ PubMed ] [ Google Scholar ] 17. Roldán E., Lisica A., Sánchez-Taltavull D., et al. Stochastic resetting in backtrack recovery by RNA polymerases. Phys. Rev. E. 2016;93:062411. doi: 10.1103/PhysRevE.93.062411. [ DOI ] [ PubMed ] [ Google Scholar ] 18. Xu P.B., Deng W.H. Lévy walk with multiple internal states. J. Stat. Phys. 2018;173:1598–1613. [ Google Scholar ] 19. Xu P.B., Deng W.H. Fractional compound Poisson processes with multiple internal states. Math. Model. Nat. Phenom. 2018;13(10) [ Google Scholar ] 20. Wang X.D., Chen Y., Deng W.H. Aging two-state process with Lévy walk and Brownian motion. Phys. Rev. E. 2019;100:012136. doi: 10.1103/PhysRevE.100.012136. [ DOI ] [ PubMed ] [ Google Scholar ] 21. Wang X.D., Chen Y., Deng W.H. Strong anomalous diffusion in two-state process with Lévy walk and Brownian motion. Phys. Rev. Res. 2020;2:013102. doi: 10.1103/PhysRevE.100.012136. [ DOI ] [ PubMed ] [ Google Scholar ] 22. Zhou T., Xu P.B., Deng W.H. Lévy walk dynamics in non-static media. J. Phys. A Math. Theor. 2022;55:025001. doi: 10.1007/s10955-022-02904-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Fedotov S., Iomin A. Migration and proliferation dichotomy in tumor-cell invasion. Phys. Rev. Lett. 2007;98 doi: 10.1103/PhysRevLett.98.118101. [ DOI ] [ PubMed ] [ Google Scholar ] 24. Masoliver J., Montero M. Anomalous diffusion under stochastic resettings: A general approach. Phys. Rev. E. 2019;100:042103. doi: 10.1103/PhysRevE.100.042103. [ DOI ] [ PubMed ] [ Google Scholar ] 25. Evans M.R., Majumdar S.N., Schehr G. Stochastic resetting and applications. J. Phys. A Math. Theor. 2020;53(19) [ Google Scholar ] 26. Lemons D.S., Gythiel A. Paul Langevin’s 1908 paper “On the Theory of Brownian Motion” [“Sur la théorie du mouvement brownien,” C. R. Acad. Sci. (Paris) 146, 530–533 (1908)] Am. J. Phys. 1997;65(11):1079–1081. [ Google Scholar ] 27. Kubo R. The fluctuation-dissipation theorem. Rep. Prog. Phys. 1966;29(1):255. [ Google Scholar ] 28. Viñales A.D., Despósito M.A. Anomalous diffusion induced by a Mittag-Leffler correlated noise. Phys. Rev. E. 2007;75:042102. doi: 10.1103/PhysRevE.75.042102. [ DOI ] [ PubMed ] [ Google Scholar ] 29. Viñales A.D., Wang K.G., Despósito M.A. Anomalous diffusive behavior of a harmonic oscillator driven by a Mittag-Leffler noise. Phys. Rev. E. 2009;80:011101. doi: 10.1103/PhysRevE.80.011101. [ DOI ] [ PubMed ] [ Google Scholar ] 30. Sandev T., Tomovski Ž., L.A. Dubbeldam J. Generalized langevin equation with a three parameter Mittag-Leffler noise. Phys. A. 2011;390(21):3627–3636. [ Google Scholar ] 31. Sandev T. Generalized langevin equation and the prabhakar derivative. Mathematics. 2017;5(4):66. [ Google Scholar ] 32. Liemert A., Sandev T., Kantz H. Generalized langevin equation with tempered memory kernel. Phys. A. 2017;466:356–369. [ Google Scholar ] 33. Helbing D., Molnár P. Social force model for pedestrian dynamics. Phys. Rev. E. 1995;51:4282–4286. doi: 10.1103/physreve.51.4282. [ DOI ] [ PubMed ] [ Google Scholar ] 34. Zhou T., Wang H., Deng W.H. Feynman-Kac equation for Brownian non-Gaussian polymer diffusion. J. Phys. A Math. Theor. 2024;57(28) [ Google Scholar ] 35. Fogedby H.C. Langevin equations for continuous time Lévy flights. Phys. Rev. E. 1994;50:1657–1660. doi: 10.1103/physreve.50.1657. [ DOI ] [ PubMed ] [ Google Scholar ] 36. Wang X.D., Chen Y., Deng W.H. Lévy-walk-like Langevin dynamics. New J. Phys. 2019;21(1):013024. [ Google Scholar ] 37. Chen Y., Wang X.D., Deng W.H. Langevin dynamics for a Lévy walk with memory. Phys. Rev. E. 2019;99:012135. doi: 10.1103/PhysRevE.99.012135. [ DOI ] [ PubMed ] [ Google Scholar ] 38. Chen Y., Wang X.D., Deng W.H. Langevin picture of Lévy walk in a constant force field. Phys. Rev. E. 2019;100:062141. doi: 10.1103/PhysRevE.100.062141. [ DOI ] [ PubMed ] [ Google Scholar ] 39. Chen Y., Deng W.H. Lévy-walk-like Langevin dynamics affected by a time-dependent force. Phys. Rev. E. 2021;103:012136. doi: 10.1103/PhysRevE.103.012136. [ DOI ] [ PubMed ] [ Google Scholar ] 40. Weron A., Magdziarz M. Anomalous diffusion and semimartingales. Europhys. Lett. 2009;86(6):60010. [ Google Scholar ] 41. Chen Z.Q., Deng W.H., Xu P.B. Feynman-Kac transform for anomalous processes. SIAM J. Math. Anal. 2021;53(5):6017–6047. [ Google Scholar ] 42. Magdziarz M. Path properties of subdiffusion-a martingale approach. Stoch. Models. 2010;26(2):256–271. [ Google Scholar ] 43. Wyłomańska A., Kumar A., Połoczański R., et al. Inverse Gaussian and its inverse process as the subordinators of fractional Brownian motion. Phys. Rev. E. 2016;94:042128. doi: 10.1103/PhysRevE.94.042128. [ DOI ] [ PubMed ] [ Google Scholar ] 44. Chen Y., Wang X.D., Deng W.H. Localization and ballistic diffusion for the tempered fractional Brownian-Langevin motion. J. Stat. Phys. 2017;169:18–37. [ Google Scholar ] 45. Chen Y., Wang X.D., Deng W.H. Resonant behavior of the generalized Langevin system with tempered Mittag-Leffler memory kernel. J. Phys. A Math. Theor. 2018;51(18) [ Google Scholar ] 46. Chen Y., Wang X.D., Deng W.H. Tempered fractional Langevin-Brownian motion with inverse β -stable subordinator. J. Phys. A Math. Theor. 2018;51(49) [ Google Scholar ] 47. Metzler R., Klafter J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 2000;339:1–77. [ Google Scholar ] 48. Barkai E., Metzler R., Klafter J. From continuous time random walks to the fractional Fokker-Planck equation. Phys. Rev. E. 2000;61:132–138. doi: 10.1103/physreve.61.132. [ DOI ] [ PubMed ] [ Google Scholar ] 49. Deng W.H., Li B.Y., Tian W.Y., et al. Boundary problems for the fractional and tempered fractional operators. Multiscale Model. Sim. 2018;16(1):125–149. [ Google Scholar ] 50. Deng W.H., Wang X.D., Zhang P.W. Anisotropic nonlocal diffusion operators for normal and anomalous dynamics. Multiscale Model. Sim. 2020;18(1):415–443. [ Google Scholar ] 51. Wu X.C., Deng W.H., Barkai E. Tempered fractional Feynman-Kac equation: Theory and examples. Phys. Rev. E. 2016;93:032151. doi: 10.1103/PhysRevE.93.032151. [ DOI ] [ PubMed ] [ Google Scholar ] 52. Wang X.D., Chen Y., Deng W.H. Feynman-Kac equation revisited. Phys. Rev. E. 2018;98:052114. [ Google Scholar ] 53. Deng W.H., Wu X.C., Wang W.L. Mean exit time and escape probability for the anomalous processes with the tempered power-law waiting times. Europhys. Lett. 2017;117(1):10009. [ Google Scholar ] 54. Chen M.H., Deng W.H. Fourth order accurate scheme for the space fractional diffusion equations. SIAM J. Numer. Anal. 2014;52(3):1418–1438. [ Google Scholar ] 55. Tian W.Y., H Z., Deng W.H. A class of second order difference approximations for solving space fractional diffusion equations. Math. Comp. 2015;84(294):1703–1727. [ Google Scholar ] 56. Zhang Z.J., Deng W.H., Karniadakis G.E. A Riesz basis Galerkin method for the tempered fractional laplacian. SIAM J. Numer. Anal. 2018;56(5):3010–3039. [ Google Scholar ] 57. Deng W.H., Li B.Y., Qian Z., et al. Time discretization of a tempered fractional Feynman-Kac equation with measure data. SIAM J. Numer. Anal. 2018;56(6):3249–3275. [ Google Scholar ] 58. Nie D.X., Deng W.H. A unified convergence analysis for the fractional diffusion equation driven by fractional Gaussian noise with Hurst index h ∈ ( 0 , 1 ) SIAM J. Numer. Anal. 2022;60(3):1548–1573. [ Google Scholar ] 59. E W., Ma C., Wu L., et al. Towards a mathematical understanding of neural network-based machine learning: What we know and what we don’t. CSIAM T. Appl. Math. 2020;1(4):561–615. [ Google Scholar ] 60. W.E, Yu B. The deep Ritz method: A deep learning-based numerical algorithm for solving variational problems. Commun. Math. Stat. 2018;6(1):1–12. [ Google Scholar ] 61. Sirignano J., Spiliopoulos K. DGM: A deep learning algorithm for solving partial differential equations. J. Comput. Phys. 2018;375:1339–1364. [ Google Scholar ] 62. Raissi M., Perdikaris P., Karniadakis G.E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys. 2019;378:686–707. [ Google Scholar ] 63. Zang Y., Bao G., Ye X., et al. Weak adversarial networks for high-dimensional partial differential equations. J. Comput. Phys. 2020;411 [ Google Scholar ] 64. Han J., Jentzen A., E W. Solving high-dimensional partial differential equations using deep learning. Proc. Natl. Acad. Sci. 2018;115(34):8505–8510. doi: 10.1073/pnas.1718942115. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 65. E W., Han J., Jentzen A. Algorithms for solving high dimensional PDEs: From nonlinear Monte Carlo to machine learning. Nonlinearity. 2021;35(1):278–310. [ Google Scholar ] 66. Han J., Long J. Convergence of the deep BSDE method for coupled FBSDEs. Probab. Uncertain. Quant. Risk. 2020;5(5) [ Google Scholar ] 67. Gao C., Gao S., Hu R., et al. Convergence of the backward deep BSDE method with applications to optimal stopping problems. SIAM J. Financ. Math. 2023;14(4):1290–1303. [ Google Scholar ] 68. Wang H., Deng W.H. Solving bivariate kinetic equations for polymer diffusion using deep learning. J. Mach. Learn. 2024;3(2):215–244. [ Google Scholar ] 69. Kenkre V.M., Montroll E.W., Shlesinger M.F. Generalized master equations for continuous-time random walks. J. Stat. Phys. 1973;9 [ Google Scholar ] 70. Scher H., Lax M. Continuous time random walk model of hopping transport: Application to impurity conduction. J. Non-Cryst. Solids. 1972;8–10:497–504. [ Google Scholar ] 71. Wang B., Anthony S.M., Bae S.C., et al. Anomalous yet Brownian. Proc. Natl. Acad. Sci. 2009;106(36):15160–15164. doi: 10.1073/pnas.0903554106. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 72. Chechkin A.V., Seno F., Metzler R., et al. Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities. Phys. Rev. X. 2017;7:021002. [ Google Scholar ] 73. Montanari A., Zecchina R. Optimizing searches via rare events. Phys. Rev. Lett. 2002;88 doi: 10.1103/PhysRevLett.88.178701. [ DOI ] [ PubMed ] [ Google Scholar ] 74. Bénichou O., Coppey M., Moreau M., et al. Optimal search strategies for hidden targets. Phys. Rev. Lett. 2005;94 doi: 10.1103/PhysRevLett.94.198101. [ DOI ] [ PubMed ] [ Google Scholar ] 75. Bartumeus F., Catalan J. Optimal search behavior and classic foraging theory. J. Phys. A Math. Theor. 2009;42(43) [ Google Scholar ] 76. Bell W.J. Springer Dordrecht; 2012. Searching Behaviour. [ DOI ] [ Google Scholar ] 77. Luby M., Sinclair A., Zuckerman D. Optimal speedup of Las Vegas algorithms. Inform. Process. Lett. 1993;47(4):173–180. [ Google Scholar ] 78. Janson S., Peres Y. Hitting times for random walks with restarts. SIAM J. Discret. Math. 2012;26(2):537–547. [ Google Scholar ] 79. Levikson B. The age distribution of Markov processes. J. Appl. Probab. 1977;14(3):492–506. [ Google Scholar ] 80. Brockwell P.J. The extinction time of a birth, death and catastrophe process and of a related diffusion model. Adv. Appl. Probab. 1985;17(1):42–52. [ Google Scholar ] 81. Kyriakidis E.G. Stationary probabilities for a simple immigration-birth-death process under the influence of total catastrophes. Stat. Probab. Lett. 1994;20(3):239–240. [ Google Scholar ] 82. Landman K., Pettet G., Newgreen D. Mathematical models of cell colonization of uniformly growing domains. B. Math. Biol. 2003;65(2):235–262. doi: 10.1016/S0092-8240(02)00098-8. [ DOI ] [ PubMed ] [ Google Scholar ] 83. Compte A. Stochastic foundations of fractional dynamics. Phys. Rev. E. 1996;53:4191–4193. doi: 10.1103/physreve.53.4191. [ DOI ] [ PubMed ] [ Google Scholar ] 84. Metzler R., Klafter J. The restaurant at the end of the random walk: Recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A Math. Gen. 2004;37(31) [ Google Scholar ] 85. Luchko Y., Gorenflo R. Scale-invariant solutions of a partial differential equation of fractional order. Fract. Calc. Appl. Anal. 1998;1 [ Google Scholar ] 86. Cairoli A., Baule A. Anomalous processes with general waiting times: Functionals and multipoint structure. Phys. Rev. Lett. 2015;115 doi: 10.1103/PhysRevLett.115.110601. [ DOI ] [ PubMed ] [ Google Scholar ] 87. Sokolov I.M., Klafter J. From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion. Chaos. 2005;15(2):026103. doi: 10.1063/1.1860472. [ DOI ] [ PubMed ] [ Google Scholar ] 88. Gorenflo R., Luchko Y., Yamamoto M. Time-fractional diffusion equation in the fractional sobolev spaces. Fract. Calc. Appl. Anal. 2015;18:799–820. [ Google Scholar ] 89. Zhou T., Trajanovski P., Xu P., et al. Generalized diffusion and random search processes. J. Stat. Mech. 2022;2022(9):093201. [ Google Scholar ] 90. Cairoli A., Baule A. Feynman-Kac equation for anomalous processes with space- and time-dependent forces. J. Phys. A Math. Theor. 2017;50(16) [ Google Scholar ] 91. Risken H. Springer Berlin, Heidelberg; 1996. The Fokker-Planck Equation. [ DOI ] [ Google Scholar ] 92. Lubich C. Discretized fractional calculus. SIAM J. Math. Anal. 1986;17(3):704–719. [ Google Scholar ] 93. Jin B., Li B., Zhou Z. Numerical analysis of nonlinear subdiffusion equations. SIAM J. Numer. Anal. 2018;56(1):1–23. [ Google Scholar ] 94. Jin B., Lazarov R., Zhou Z. Numerical methods for time-fractional evolution equations with nonsmooth data: A concise overview. Comput. Methods Appl. Mech. Eng. 2019;346:332–358. [ Google Scholar ] 95. Sun J., Nie D.X., Deng W.H. Fast algorithms for convolution quadrature of Riemann-Liouville fractional derivative. Appl. Numer. Math. 2019;145:384–410. [ Google Scholar ] 96. Podlubny I. Elsevier; 1998. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, Methods of Their Solution and Some of Their Applications. [ Google Scholar ] 97. Sabzikar F., Meerschaert M.M., Chen J. Tempered fractional calculus. J. Comput. Phys. 2015;293:14–28. doi: 10.1016/j.jcp.2014.04.024. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 98. Acosta G., Borthagaray J.P. A fractional Laplace equation: Regularity of solutions and finite element approximations. SIAM J. Numer. Anal. 2017;55(2):472–495. [ Google Scholar ] 99. Bonito A., Borthagaray J.P., Nochetto R.H., et al. Numerical methods for fractional diffusion. Comput. Vis. Sci. 2018;19(5):19–46. [ Google Scholar ] 100. Jia R.-Q. Spline wavelets on the interval with homogeneous boundary conditions. Adv. Comput. Math. 2009;30(2):177–200. [ Google Scholar ] 101. Song R., Vondraček Z. Potential theory of subordinate killed Brownian motion in a domain. Probab. Theory Relat. Fields. 2003;125:578–592. [ Google Scholar ] 102. Bernadotte A., Mikhelson V.M., Spivak I.M. Markers of cellular senescence. Telomere shortening as a marker of cellular senescence. Aging. 2016;8(1):3–11. doi: 10.18632/aging.100871. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 103. Lingner J., Hughes T.R., Shevchenko A., et al. Reverse transcriptase motifs in the catalytic subunit of telomerase. Science. 1997;276(5312):561–567. doi: 10.1126/science.276.5312.561. [ DOI ] [ PubMed ] [ Google Scholar ] 104. Han P.P., Zhou Y., Deng W.H. Modeling telomere shortening process. Quant. Biol. 2025;13(1) doi: 10.1002/qub2.74. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 105. Alder J.K., Hanumanthu V.S., Strong M.A., et al. Diagnostic utility of telomere length testing in a hospital-based setting. Proc. Natl. Acad. Sci. 2018;115(10):E2358–E2365. doi: 10.1073/pnas.1720427115. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 106. Golding I., Cox E.C. Physical nature of bacterial cytoplasm. Phys. Rev. Lett. 2006;96:098102. doi: 10.1103/PhysRevLett.96.098102. [ DOI ] [ PubMed ] [ Google Scholar ] 107. Molina-Garcia D., Sandev T., Safdari H., et al. Crossover from anomalous to normal diffusion: Truncated power-law noise correlations and applications to dynamics in lipid bilayers. New J. Phys. 2018;20(10) [ Google Scholar ] 108. Janczura J., Orzeł S., Wyłomańska A. Subordinated α -stable Ornstein-Uhlenbeck process as a tool for financial data description. Phys. A. 2011;390(23):4379–4387. [ Google Scholar ] 109. Scher H., Margolin G., Metzler R., et al. The dynamical foundation of fractal stream chemistry: The origin of extremely long retention times. Geophys. Res. Lett. 2002;29(5) [ Google Scholar ] 110. Nezhadhaghighi M.G., Rajabpour M.A., Rouhani S. First-passage-time processes and subordinated Schramm-Loewner evolution. Phys. Rev. E. 2011;84:011134. doi: 10.1103/PhysRevE.84.011134. [ DOI ] [ PubMed ] [ Google Scholar ] 111. Morgan D.O. New Science Press; 2007. The Cell Cycle: Principles of Control. [ Google Scholar ] 112. Gudimchuk N.B., McIntosh J.R. Regulation of microtubule dynamics, mechanics and function through the growing tip. Nat. Rev. Mol. Cell Biol. 2021;22:777–795. doi: 10.1038/s41580-021-00399-x. [ DOI ] [ PubMed ] [ Google Scholar ] 113. Shaevitz J.W., Abbondanzieri E.A., Landick R., et al. Backtracking by single RNA polymerase molecules observed at near-base-pair resolution. Nature. 2003;426:684–687. doi: 10.1038/nature02191. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 114. Zheng X. Springer Berlin, Heidelberg; 2009. Mechanics of Wind-Blown Sand Movements. [ DOI ] [ Google Scholar ] 115. Kok J.F., Parteli E.J.R., Michaels T.I., et al. The physics of wind-blown sand and dust. Rep. Prog. Phys. 2012;75(10) doi: 10.1088/0034-4885/75/10/106901. [ DOI ] [ PubMed ] [ Google Scholar ] 116. Brunet Y. Turbulent flow in plant canopies: Historical perspective and overview. Bound. -Lay. Meteorol. 2020;177:315–364. [ Google Scholar ] 117. de Langre E. Effects of wind on plants. Annu. Rev. Fluid Mech. 2008;40:141–168. [ Google Scholar ] 118. Price J.F. Upper ocean response to a hurricane. J. Phys. Oceanogr. 1981;11(2):153–175. [ Google Scholar ] 119. Liu L., Fei J.F., Cheng X.P., et al. Effect of wind-current interaction on ocean response during Typhoon KAEMI (2006) Sci. China Earth Sci. 2013;56:418–433. [ Google Scholar ] 120. Kwak K., Song H., Marshall J., et al. Suppressed pCO2 in the southern ocean due to the interaction between current and wind. J. Geophys. Res. Oceans. 2021;126(12) [ Google Scholar ] 121. England M.H., McGregor S., Spence P., et al. Recent intensification of wind-driven circulation in the pacific and the ongoing warming hiatus. Nat. Clim. Chang. 2014;4:222–227. [ Google Scholar ] 122. Mazoyer C., Vanneste H., Dufresne C., et al. Impact of wind-driven circulation on contaminant dispersion in a semi-enclosed bay. Estuar. Coast. Shelf S. 2020;233 [ Google Scholar ] 123. Stevens A. A stochastic cellular automaton modeling gliding and aggregation of myxobacteria. SIAM J. Appl. Math. 2000;61(1):172–182. [ Google Scholar ] Articles from Fundamental Research are provided here courtesy of The Science Foundation of China Publication Department, The National Natural Science Foundation of China ACTIONS View on publisher site PDF (955.7 KB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top