## Abstract

Single-particle tracking offers detailed information about the motion of molecules in complex environments such as those encountered in live cells, but the interpretation of experimental data is challenging. One of the most powerful tools in the characterization of random processes is the power spectral density. However, because anomalous diffusion processes in complex systems are usually not stationary, the traditional Wiener-Khinchin theorem for the analysis of power spectral densities is invalid. Here, we employ a recently developed tool named aging Wiener-Khinchin theorem to derive the power spectral density of fractional Brownian motion coexisting with a scale-free continuous time random walk, the two most typical anomalous diffusion processes. Using this analysis, we characterize the motion of voltage-gated sodium channels on the surface of hippocampal neurons. Our results show aging where the power spectral density can either increase or decrease with observation time depending on the specific parameters of both underlying processes.

## Introduction

A very large class of biological and physical systems exhibit correlations that extend across multiple time scales. This feature is also found in social networks as well as in complex systems made of interacting components like glasses. Such correlations manifest themselves as a broad spectrum of relaxation times and in the practically universal emergence of 1/*f* decay in the power spectrum, which points to self-similarity in the dynamics at different timescales^{1,2}. The effect is predominantly found at low frequencies where the contributions of each frequency *ω* = 2*π**f* to the overall power spectral density (PSD) exhibit a power law *S*(*ω*) ~ 1/*ω*^{β}, with 0 < *β* ≤ 2^{3,4,5,6,7}. To name a few diverse examples, 1/*f* spectra are observed in nanoscale devices^{8,9}, network traffic^{10}, earthquakes^{11}, heartbeat dynamics^{12}, DNA base sequences^{13}, climate^{14}, and ecology^{15}. Mandelbrot and later Bouchaud et al. suggested that the processes involved are inherently non-stationary leading to the idea that the spectrum should depend both on the frequency and the measurement time^{5,16,17}. Indeed, the very basic formula describing these ubiquitous phenomena was recently replaced with a more general one^{6}. Based on experimental data of blinking quantum dots^{18}, nanoelectronic devices^{9,19}, and fluctuations of interfaces^{7}, the basic spectrum must be described with a new formula \(S(\omega,{t}_{{{{{{{{\rm{m}}}}}}}}} ) \sim {\omega }^{-\beta }{t}_{{{{{{{{\rm{m}}}}}}}}}^{{{{{{\rm{z}}}}}}}\), where *t*_{m} is the measurement time. These developments, in turn, motivated a new theoretical framework, called aging Wiener–Khinchin theorem^{20,21,22}. This new theorem replaces the celebrated Wiener–Khinchin theorem valid for stationary processes, which is widely applicable to systems that do not exhibit 1/*f* noise^{23}.

Notwithstanding previous advances, many questions remain open. First, the aging Wiener–Khinchin theorem relates the aging power spectrum with *z* ≠ 0 to a non-stationary correlation function (soon to be discussed). However, how can one find this correlation function? As for the standard Wiener–Khinchin theorem, the correlation function is specific to the system. In the context of diffusion in cells as well as in many other complex systems, Mandelbrot’s fractional Brownian motion (fBM)^{24} and the Montroll–Weiss continuous time random walk (CTRW)^{25} are two widely investigated models of anomalous transport. While the fluctuations in fBM are stationary, the CTRW process is inherently non-stationary. However, both models, when standing alone, are usually non-sufficient to describe the transport of particles that alternate between a trapping phase (like in CTRW) and correlated motion (like in fBM), as is the case in live cells, for example due to interactions in a viscoelastic medium^{26}. The open questions begin with how to create a marriage between these models? Then, can we obtain the correlation functions and 1/*f* spectrum? Achieving these goals will show how the exponents *β* and *z* depend on the underlying processes, and will determine which of the processes is dominating the PSD. Finally, most importantly, these goals can elucidate whether the whole approach to the PSD is useful in experiments. Specifically, we demonstrate the applicability of aging Wiener–Khinchin theorem and the corresponding calculation of the correlation function with experimental recording of the power spectra of the motion of ion channels in the plasma membrane of mammalian cells.

The emergence of 1/*f* noise has triggered notable interest in biological environments both from a fundamental point of view and due to its relevance in pathologies and disease^{27}. Self-similar temporal characteristics are observed in biological systems of broadly different length scales. Recent molecular dynamics simulations in combination with previous experimental results have shown that the internal dynamics in globular proteins are self-similar and the autocorrelation function is aging over an astonishing 13 decades in time^{2,28}. These fluctuations play essential roles in cell functions that involve molecular interactions such as gene regulation. In fact, this behavior is widespread and found from the dynamics of proteins within cell membranes to the scaling behavior of heartbeat time series^{27,29}. Nevertheless, it still remains a challenge to measure how aging affects the spectrum of recorded 1/*f* noise in real systems.

Single molecule tracking in the cell environment has been used extensively to shed light on the functions and interactions of the molecules that make life possible^{30,31,32,33}. Spectral analyses are emerging as a key tool in the characterization of individual molecule trajectories in biological systems because it informs on features that are difficult to infer using other traditional statistics^{6,34,35,36,37,38}. It has been observed that among traditional statistical approaches, e.g., analyses based on the mean squared displacement, the PSD appears to be less sensitive to external noises^{39}. Following previous work, we promote a theory that shows how the most basic formula of 1/*f* noise needs modifications, namely that \(S(\omega,{t}_{{{{{{{{\rm{m}}}}}}}}} ) \sim {\omega }^{-\beta }{t}_{{{{{{{{\rm{m}}}}}}}}}^{{{{{{\rm{z}}}}}}}\) as mentioned. The question that still needs to be addressed is what the physical meaning of the new exponent *z* is, to explore cases where it is negative (corresponding to a decrease of the PSD with time and, hence, aging) and cases where it is positive (corresponding to a PSD increasing with time and, hence, rejuvenation). Further, beyond the development of the theory, it is important to show how these effects are found experimentally.

Traditionally, the PSD of a time-dependent signal is defined as the average over an infinitely large ensemble in the limit of infinite time (Supplementary Eq. 1). In practice, when analyzing either experiments or numerical simulations, one does not have access to infinite measurement time, nor to a large ensemble of trajectories, and the PSD is estimated by using the periodogram. For stationary processes, the PSD can be directly calculated from the autocorrelation function, using the relation provided by the Wiener–Khinchin theorem (Supplementary Eq. 2)^{23}. The Wiener–Khinchin theorem holds for a large class of time-invariant processes, where the concept of a time-independent limiting power spectrum is useful. One could wonder how to extend the Wiener–Khinchin theorem to non-stationary processes, but, due to the extensive variety of such processes, this general approach appears a priori to be a futile direction of research. Nonetheless, this first assessment turns out to be wrong. There exists a large class of stochastic processes describing systems that are non-stationary but scale invariant. Specifically, the autocorrelation function explicitly depends on time *t* via the expression *C*_{EA}(*τ*, *t*) = 〈*x*(*t*) *x*(*t* + *τ*)〉 ~ *t*^{γ}*ϕ*_{EA}(*τ*/*t*), where *ϕ*_{EA}(*τ*/*t*) is a scaling function. As mentioned, a new theoretical framework was developed for this very large class of scale invariant processes, the aging Wiener–Khinchin theorem^{20,21,22}. The PSD that emerges in this case is, in turn, directly related to 1/*f* noise and depends on the observation time.

Here, we address the spectral content of processes with scale free relaxation times, using both theoretical modeling and experimental validation. We show how the aging Wiener–Khinchin theorem is a useful tool and, more importantly, demonstrate how the aging exponent *z* and the spectral exponent *β* are related to the underlying processes. To reach this goal, we obtain the non-stationary correlation function of the subordinated fBM, which combines two well known approaches to anomalous diffusion. Depending on whether the process is negatively or positively correlated, we get vastly different frequency decays of the power spectrum. Thus, the aging Wiener–Khinchin theorem can be used to classify widely different classes of dynamics. Finally, by analyzing the dynamics of voltage-gated sodium channels (Nav) on the somatic membrane of hippocampal neurons, we demonstrate the usefulness of the approach, and prove that its basic principles work in the laboratory. These experimental data reveal how one can use a few long trajectories and estimate the exponents characterizing the dynamics with high precision. Our work, thus, not only validates the aging Wiener-Khinchin theorem as an emerging tool in spectral analysis, but it also unravels the meaning of the exponents describing the aging and the frequency decay.

## Results

### Aging Wiener–Khinchin theorem

In any stationary process *x*(*t*), the PSD is related to the autocorrelation function (ACF) *C*_{EA}(*τ*) = 〈*x*(*t*)*x*(*t* + *τ*)〉 via the fundamental Wiener–Khinchin theorem (Supplementary Eq. 2). Throughout the manuscript we employ the subscripts EA and TA to denote ensemble averages and time averages, respectively. However, diffusive processes are intrinsically non-stationary and thus the Wiener–Khinchin theorem is invalid. In recent years, power spectrum theory has been expanded with a tool called the aging Wiener–Khinchin theorem^{20,21,22}. This theorem covers a broad class of non-stationary processes that possess an autocorrelation function with the long-time asymptotic *C*_{EA}(*t*, *τ*) = 〈*x*(*t*)*x*(*t* + *τ*)〉 ~ *t*^{γ}*ϕ*_{EA}(*τ*/*t*). Such correlation functions are common^{22,40,41} and they are called scale invariant. An alternative analysis of the autocorrelation function is performed in terms of its time average *C*_{TA} of individual trajectories, where

with *t*_{m} being the measurement time. For ergodic processes, *C*_{TA} converges to *C*_{EA} in the long time limit. However, when the process is not ergodic, such as a scale-free CTRW, *C*_{TA} of individual trajectories remain random variables even in the long time limit^{42,43}. Thus, one analyzes the ensemble-average of the TA-ACF, 〈*C*_{TA}(*t*_{m}, *τ*)〉. Further, ergodicity breaking leads to a difference in the two averages, 〈*C*_{TA}(*t*_{m} = *t*, *τ*)〉 ≠ *C*_{EA}(*t*, *τ*). Each of these formalisms (ensemble vs. time averages) has its own advantages and disadvantages. Nevertheless, when the number of trajectories is small and the measurement time is long, the time averages lead to better statistics and it is, thus, the more commonly used method in single-particle tracking. When *C*_{EA}(*t*, *τ*) = *t*^{γ}*ϕ*_{EA}(*τ*/*t*), the time-average ACF has also the scaling form \(\langle {C}_{{{{{{{{\rm{TA}}}}}}}}}({t}_{{{{{{{{\rm{m}}}}}}}}},\tau )\rangle ={t}_{{{{{{{{\rm{m}}}}}}}}}^{\gamma }{\phi }_{{{{{{{{\rm{TA}}}}}}}}}(\tau /{t}_{{{{{{{{\rm{m}}}}}}}}})\)^{20}. The scaling function *ϕ*_{TA}(*τ*/*t*_{m}) is directly related to the ensemble average via the relation

where *y* = *τ*/*t*_{m}, which implies 0 ≤ *y* ≤ 1.

For a measurement time *t*_{m} the power spectrum can be only obtained for the discrete set of frequencies *ω*_{k}*t*_{m} = 2*π**k* with *k* being a non-negative integer. That is, the frequencies can be resolved down to Δ*ω* = 2*π*/*t*_{m}, which decays to zero in the limit of large measurement time *t*_{m}. The aging Wiener–Khinchin theorem relates the average power spectrum for this set of frequencies to the time-averaged autocorrelation function^{20,22},

A relation between the PSD and the ensemble-averaged correlation function also exists, but we will employ the relation to the time average because of its more common use in single-particle tracking experiments.

### The model for subordinated random walks

A useful way to model the diffusive transport in live cells is via the combination of two stochastic processes: the CTRW and fBM. On one hand, the CTRW constitutes the quintessential diffusion process with heavy-tailed immobilization times and has been extensively used to describe transport in disordered environments^{44,45}, protein dynamics in mammalian cells^{30,31,32,46}, and even to model financial markets^{47}. On the other hand processes with correlated increments such as fBM or diffusion in fractal environments have been often observed to lead to anomalous transport with memory effects^{48,49,50}. fBM is the only Gaussian self-similar process with stationary increments, of which Brownian motion constitutes a special case. Technically the combination of these widely observed models is made possible with a subordination technique^{32,51,52}. In a subordination scheme, the steps of a random walk take place at operational times *t*_{n} defined by a directing stochastic process. For example, antipersistent motions accompanied by heavy-tailed immobilization times, have been observed in live cells in the motion of ion channels^{53}, insulin granules^{54}, membrane receptors^{55}, and nanosized objects in the cytoplasm^{56}, as well as for tracer particles in actin networks in vitro^{57}. Subordinated processes constitute one of the most general classes of random walks and are widespread beyond the dynamics in the cell^{4,7,19,58}. This scheme allows to evaluate processes with short-range or long-range memory and non-stationarity, leading to complex aging properties.

We consider a fBM-like process at discrete times, *n* = 0, 1, 2, 3,…, with Hurst exponent *H*, such that its autocorrelation function at the discrete times *n* is given by^{24}

where the coefficient Δ*x* is a scaling parameter with units of length. We place the process defined by Eq. (4) under the operational time of a CTRW, so that the particle is immobilized during sojourn times with a heavy-tailed distribution. Such immobilizations arise, for example, from energetic disorder where a particle has random waiting times at each trapping site^{25,29,59,60}.

The operational times are defined by a random process {*t*_{n}} with non-negative independent increments *τ*_{n} = *t*_{n} − *t*_{n−1}. The time increments *τ*_{n} between renewals are, in the long time limit, asymptotically distributed according to a probability density function^{61}

where 0 < *α* < 1, *t*_{0} is a constant with units of time, and Γ(⋅) is the gamma function. At time *t*, the position of the particle is *x*(*t*) = *x*_{n} where *n* is the random number of renewals in the interval (0, *t*). Given *n*, the position *x*_{n} is determined by the discrete fBM process defined by Eq. 4. Three representative trajectories of such a process are shown in Fig. 1. The ensemble-averaged autocorrelation function of *x*(*t*) is then

where \({\mathbb{E}}[g(x)]=\langle g(x)\rangle\) represents the expected value of *g*(*x*) and \({\mathbb{E}}[g(x)| y]\) is the conditional expected value of *g*(*x*) given *y*. In particular, the last term indicates the iterated expectation of *x*(*t*)*x*(*t* + *τ*), given that *n* steps have taken place up to time *t* and *n* + Δ*n* steps have taken place up to time *t* + *τ*. Further, we define *χ*_{n,Δn}(*t*, *τ*) as the joint probability of taking *n* steps up to time *t* and Δ*n* steps in the interval (*t*, *t* + *τ*). Combining Eqs. (6) and (4), we obtain

Once the ensemble-averaged autocorrelation function is found, we can obtain the time-averaged *C*_{TA}(*t*_{m}, *τ*) via Eq. (2) and, subsequently, the PSD using the aging Wiener–Khinchin theorem (Eq. 3).

### Continuous time random walk (2*H* = 1)

The fBM reverts to Brownian motion when 2*H* = 1 and, thus, the process becomes a traditional CTRW^{25,62}. The ensemble-averaged autocorrelation function in Eq. (7) becomes (see Supplementary Note 3)

which, given the memoryless property of Brownian motion, boils down to the ensemble-averaged autocorrelation function being independent of lag time *τ* and equal to the mean squared displacement (MSD), *C*_{EA}(*t*, *τ*) = 2Δ*x*^{2}〈*n*(*t*)〉 = 〈*x*^{2}(*t*)〉. The MSD solution for the CTRW is 〈*x*^{2}(*t*)〉 ~ *t*^{α}, that is, it exhibits subdiffusion with anomalous exponent *α*^{61}.

The ensemble-averaged autocorrelation function in Eq. (8), for 2*H* = 1, implies that *C*_{EA} = *t*^{α}*ϕ*_{EA} with *ϕ*_{EA} being a constant. The time-averaged autocorrelation function is \(\langle {C}_{{{{{{{{\rm{TA}}}}}}}}}\rangle ={t}_{{{{{{{{\rm{m}}}}}}}}}^{\alpha }{\phi }_{{{{{{{{\rm{TA}}}}}}}}}(\tau /{t}_{{{{{{{{\rm{m}}}}}}}}})\) and we find (Supplementary Eq. 17)

Next, we use the time-averaged autocorrelation function (Eq. 9) in conjunction with the aging Wiener–Khinchin theorem to obtain the PSD of the CTRW. We find the exact solution of the sample power spectral density by solving the integral in Eq. (3). The PSD (Supplementary Eq. 18) is a function of both frequency *ω* and realization time *t*_{m}. Expanding the PSD for *ω**t*_{m} ≫ 1, it is found that the leading term scales in frequency as *ω*^{−2} and in time as \({t}_{{{{{{{{\rm{m}}}}}}}}}^{-(1-\alpha )}\),

which is related to the MSD via the relation

This is a useful relation that connects the fluctuations in the trajectory (the PSD) to transport properties (the MSD) for the CTRW. Importantly, the MSD is proportional to the mean number of renewals, thus Eq. (11) provides a connection between the PSD and the number of renewals. While Eq. (11) applies to the CTRW, we will see later that it is not universal for the scale free processes under study.

Figure 2 shows a comparison of these analytical results to numerical simulations of 10,000 realizations with *α* = 0.7. The MSD exhibits a power law, 〈*x*^{2}(*t*)〉 ~ *t*^{α} (Fig. 2a). The power spectral density is presented in Fig. 2b for five different measurement times from *t*_{m} = 2^{8} to 2^{16} and shows good agreement with the power law asymptotic *ω*^{−2}. As shown in Supplementary Note 3, using hypergeometric functions we can get the exact PSD; however, the power law asymptotics show highly accurate results. The spectra also exhibit aging with an amplitude that scales as \({t}_{{{{{{{{\rm{m}}}}}}}}}^{-(1-\alpha )}\) (Fig. 2c). Intuitively, as the measurement time increases, we encounter longer stagnation periods and, hence, the PSD decays with measurement time. Physically, this effect is due to the broadly distributed trapping times in the system.

### Subordinated process involving fBM (0 < *H* < 1)

We now deal with subordinated random walks where the increments exhibit correlations. When 2*H* ≠ 1, the process has positively correlated increments for *H* > 0.5 and negatively correlated increments for *H* < 0.5. The autocorrelation function *C*_{EA} in Eq. (7) is

where Δ*n*(*τ*; *t*) is the number of steps between the aged time *t* and *t* + *τ*. Using renewal theory and the power law waiting time distribution in Eq. (5), the terms in Eq. (12) can be expressed via hypergeometric functions (Supplementary Eqs. 21 and 22). The ensemble-averaged autocorrelation function (Supplementary Eq. 25) has the form *C*_{EA}(*t*, *τ*) = *t*^{γ}*ϕ*_{EA}(*τ*/*t*), which implies the time-averaged autocorrelation function is of the form \(\langle {C}_{{{{{{{{\rm{TA}}}}}}}}}({t}_{{{{{{{{\rm{m}}}}}}}}},\tau )\rangle ={t}_{{{{{{{{\rm{m}}}}}}}}}^{\gamma }{\phi }_{{{{{{{{\rm{TA}}}}}}}}}(\tau /{t}_{{{{{{{{\rm{m}}}}}}}}})\)^{22}. Following Eq. (2), we find the scaling function *ϕ*_{TA}(*τ*/*t*_{m}). The exact analytical results for the time-averaged ACF (Supplementary Eq. 27) were compared to numerical simulations. The simulations are observed to agree with analytical results for both *H* < 0.5 and *H* > 0.5 in Supplementary Fig. (1a, b), respectively.

The calculation of the PSD with the correlation function involves two steps. Our approach uses the scale invariant correlation function, which was tested versus numerical results, and the aging Wiener–Khinchin theorem, Eq. (3). The calculation essentially leads to PSDs that are expressed in terms of hypergeometric functions (Supplementary Eq. 33) and can be simplified. The idea is to use the large frequency limit to obtain approximate results of the aging 1/*f* noise type. These work well, as we show later in the figures. By expanding the PSD in the limit *ω**t*_{m} ≫ 1 and noting that the spectrum is evaluated at frequencies *ω**t*_{m} = 2*π**k*, we obtain the leading term, which depends on the specific values of *α* and *H*. In the case that the increments are anticorrelated, i.e., *H* < 0.5,

where *c* is a constant defined explicitly in Supplementary Eq. (38). An example of this antipersistent case is shown for numerical simulations with *α* = 0.4 and *H* = 0.3 in Fig. 3a. The scaling of the PSD both in *t*_{m} and *ω* agrees with Eq. (13).

When the increments of the random walk are positively correlated (i.e, *H* > 1/2), the leading term is

with *D* being a generalized diffusion coefficient (Supplementary Eq. 26). This PSD is related to the mean square displacement in a similar way as the CTRW, via the relation

which is similar to Eq. (11), albeit with a factor 1/2. When *H* > 1/2, the PSD decreases with observation time for small *α* and *H*, namely when 2*α**H* < 1. Otherwise (shaded regime III in Fig. 3d), the PSD increases with observation time. Figure 3b shows the power spectra for numerical simulations where the underlying fBM is superdiffusive with *H* = 0.7 and *α* = 0.4, which falls in the regime that 〈*S*(*ω*, *t*_{m})〉 decays with *t*_{m} (regime II in Fig. 3d). Figure 3c shows simulations with *H* = 0.75 and *α* = 0.8 where 〈*S*(*ω*, *t*_{m})〉 indeed is observed to increase with *t*_{m}. In this regime of increasing *S*, the convergence to Eq. (14) is very slow and appears to converge only for realization times *t*_{m} > 10^{5}. The increase of 〈*S*(*ω*, *t*_{m})〉 with time is directly related to the persistent property of the fBM^{36}.

We now focus on two important limits of our results, namely, the limits 2*H* → 1 and *α* → 1. In the first one, the process reverts to the traditional CTRW, for which the result is given by Eq. (10). Here, the two leading terms in the exact result for the PSD (Supplementary Eq. 33) converge to the same exponent yielding the simple asymptotic approximation 〈*S*_{2H=1}(*ω*, *t*_{m})〉 ≈ 2(*D* + *c*)/*ω*^{2}. The agreement with Eq. (10) serves as a basic test to evaluate the results. The second limit (*α* → 1) is expected to converge to the known results for the standard fBM. In this limit, the PSD becomes (i) 〈*S*(*ω*, *t*_{m})〉 ~ 1/*ω*^{1+2H} when 2*H* < 1 and (ii) \(\langle S(\omega ,{t}_{{{{{{{{\rm{m}}}}}}}}})\rangle \sim {t}_{{{{{{{{\rm{m}}}}}}}}}^{2H-1}/{\omega }^{2}\) when 2*H* > 1. These expressions are in agreement with the known formulas for subdiffusive and superdiffusive fBM, respectively (see e.g.,^{36}).

### Experimental results

The derivation of the PSD of subordinated random walks enables us to characterize the motion of membrane proteins that typically interact with heterogeneous partners. These trajectories are obtained using single molecule tracking of labeled proteins in living cells. An example of a transmembrane protein that exhibits heterogeneous interactions is the voltage gated sodium channel Nav1.6. It was previously found that in the somatic plasma membrane of hippocampal neurons, Nav1.6 channels are transiently confined into cell surface nanodomains^{63}. Because these nanodomains are only of the order of 100 nm in size, we can neglect the motion within an individual domain without altering the long time statistics of the process. Further, it was reported that the motion of these channels displays ergodicity breaking due to their transient confinement^{64}. These effects lead to the idea of trapping and the CTRW type of dynamics. Thus, we model the confinement (immobilization) times using Eq. (5). An important property of heavy-tailed renewal processes is that they depend on the time that lapsed since the system started^{65}. In the case of single molecule tracking of Nav channels, measurements start when the channel is delivered to the plasma membrane and, thus, the time *t* = 0 is well-defined. Besides transient immobilizations, Nav1.6 also show antipersistent fBM-like motion, leading to a non-linear time-averaged MSD. Here, we evaluate 87 Nav1.6 trajectories of 256 data points each, with a sampling time Δ*t* = 50 ms.

Before digging into the PSD analysis of Nav channels, we consider their mean square displacement, which is a familiar statistical tool that helps us understand some basic properties of their motion. Furthermore, we can evaluate the validity of our model for the motion of membrane proteins by analyzing the relations between the exponents that characterize the mean squared displacement and the power spectrum. Figure 4a shows the ensemble-averaged MSD (EA-MSD, 〈*x*^{2}(*t*)〉) together with its 95% confidence interval and the ensemble-average of the time-averaged MSD (EA-TA-MSD) for three different observation times, *t*_{m} = 64Δ*t*, 128Δ*t*, and 256Δ*t*. The EA-TA-MSD is defined in its usual way,

where, using the same notation as in the autocorrelation function, *τ* denotes the lag time. The difference between the EA-TA-MSD and the EA-MSD (Fig. 4a) is a direct indication of ergodicity breaking in the motion of Nav channels^{30,64}. In the context of our model, the ergodic hypothesis breaks down since *α* < 1. In theory, it should be possible to use the ensemble-averaged MSD to extract information about the exponents that characterize the motion. However, when the number of trajectories is not very large (as is usually the case in live cell experiments), the estimation of exponents from this metric is very poor due to statistical errors. This effect can be directly seen in the confidence interval of the MSD in Fig. 4a. Thus, we propose here to employ in addition to the TA-MSD a robust metric such as the PSD.

The EA-TA-MSD of the subordinated process scales as^{66}

We have measured both the EA-TA-MSD (Fig. 4a) and the PSD (Fig. 4b), with different observation times *t*_{m}. From the MSD, using Eq. (17), we extract exponents *α* = 0.54 ± 0.02 and *H* = 0.32 ± 0.08. Remarkably, this is nearly identical to the estimation based on the PSD, where, using Eq. (13), we obtain *α* = 0.50 ± 0.02 and *H* = 0.25 ± 0.11. The agreement is not a coincidence and it indicates that the underlying model of a subordinated process is consistent with two independent measurements. In other words, we can use one set of measurements (e.g., PSD) to predict the exponents of the other (e.g., MSD) and show that the selected model works. From a single set of data we cannot make this conclusion. Namely, if we record *β* and *z*, we can easily estimate the exponents *α* and *H*, but that, as a stand alone, is not informative, since the number of fitting parameters (two) is the same as the number of linear equations given in the relations between the exponents (*β*, *z*) and the exponents (*α*, *H*). Hence, extraction of these exponents with an additional measurement is required to find a consistent theory, beyond merely fitting parameters.

A key aspect of these measurements is that the PSD is obtained for different observation times. By increasing the measurement time, we indeed observe the aging power spectrum, an effect that could have been missed, as the natural tendency in experiments is simply to use the longest available trajectories. The PSD decays with observation time, i.e., *z* < 0, as predicted for a process with Hurst exponent *H* < 1/2 (see Eq. 13). The PSD amplitude as a function of measurement time *t*_{m} is shown in the inset of Fig. 4b, indicating *z* = − 0.50 ± 0.02. This spectral analysis in combination with the MSD confirms the predictions stating that the motion of Nav channels is a subordinated process and lets us obtain accurate estimates of the waiting time distribution and the Hurst exponent. While the goal of this work pertained to the dynamics of proteins, it is directly applicable to any process where a correlated random walk coexists with a non-ergodic CTRW.

## Discussion

We characterize subordinated random walks via two exponents, the Hurst exponent *H* and the exponent that describes the heavy-tailed waiting time distribution *α*. The Hurst exponent governs the correlations between increments and the memory effects of the random walk, while the exponent *α* is responsible for the long waiting times. We observe that the PSD is found to be accurately described by the formula \(S(\omega ,{t}_{{{{{{{{\rm{m}}}}}}}}}) \sim {\omega }^{-\beta }{t}_{{{{{{{{\rm{m}}}}}}}}}^{z}\), where the exponents *β* and *z* are uniquely defined by *H* and *α* (see Eqs. 13 and 14).

Our results can be divided into two large classes depending on whether the increments are positively or negatively correlated, that is *H* > 1/2 or *H* < 1/2. The case *H* < 1/2 is associated with the tracer’s interactions with a viscoelastic medium which lead to subdiffusion, while *H* > 1/2 is associated with persistent walks that can lead to superdiffusion, which is, in turn, related to active transport. Let us discuss first our results for antipersistent random walks (*H* < 1/2) because this is the relevant regime for the dynamics of the membrane proteins we studied. In this situation, *β* = 2 − *α* + 2*α**H*, i.e, it is influenced by both the properties of the fBM and the CTRW, and it falls in the range 1 < *β* < 2. In contrast, the exponent *z* is dictated solely by the power law trapping times of the CTRW and it shows aging, that is *z* < 0. Specifically, the PSD decays with measurement time with an exponent *z* = *α* − 1. As such, the aging process in this regime does not contain any information about the fBM. For the Nav1.6 ion channels we recorded aging power spectra with *z* = −0.50 and *β* = 1.75, from which we estimated the Hurst index *H* of the fBM-like process and the exponent *α* that characterizes the tail of the distribution of immobilization times. Then, it is possible to use the measured exponents *z* and *β* (which give *α* and *H*) to predict the exponents of the time- and ensemble-averaged MSD, and compare these predictions to the experimental data. An agreement between predicted and measured exponents would show that the model is working well without any fitting to the MSD measurements. Indeed, the experimental results with Nav channels provide a very strong validation of our hypothesis, in which these channels can be described by an antipersistent random walk (*H* < 1/2) in the presence of traps caused by interactions with heterogeneous partners at the plasma membrane (*α* < 1). The antipersistent walk is a signature of spatial heterogeneity and self similar obstructions in the membrane, while the heavy tailed waiting times are caused by trapping events, e.g. energy disorder. Our findings that subordinated fBM is the relevant model for ion channels is significant. The slow dynamics can rationalize the organization of these membrane proteins, as practically immobile, but still have some dynamics which is important for allowing interactions with cytoskeletal components and other reaction partners. We expect that our model can be used to determine diffusion-limited reaction rates.

The case of positively correlated increments *H* > 1/2 leads to much richer phenomena, and we encounter both aging (a decay of the fluctuations) and rejuvenation (increase of the fluctuations) with measurement time. In this situation, *β* is a constant *β* = 2 independent of the exponents of either the fBM or the CTRW. Note that this is the same frequency scaling as that of Brownian motion and the traditional CTRW. Nevertheless, the exponent *z*, given by Eq. (14), presents intriguing properties and, hence, it is the aging that informs about interesting physical effects. The smaller *α*, the faster the fluctuations are inhibited over time. This effect is due to the particles becoming more and more immobile, i.e., they find deeper and deeper traps the longer the time that lapses since the preparation of the setup. However, when *H* is increased, the fBM becomes more superdiffusive and, strikingly, the PSD can be observed to rejuvenate, i.e., the fluctuations become more prominent over time. Precisely, the turnover from aging to rejuvenating takes place when *H* > 1/(2*α*). Thus, one can infer the region in phase space to which the process belongs by noting whether the fluctuations increase or decay (see Fig. 3d for a full phase diagram). The special case *z* = 0, that is so often tacitly implied in the 1/*f* literature, is actually rare in subordinated diffusive processes and it only takes place when *H* = 1/(2*α*). Note, however, that normal Brownian motion takes place in the limit that *H* = 1/2 and *α* = 1, which also implies *z* = 0, namely the absence of an aging effect.

Our analytical results describe the spectral content of a wide class of non-stationary processes with scale invariant correlation functions. The derivations are obtained using the aging Wiener–Khinchin theorem and we demonstrate the applicability of this theory with experimental trajectories of molecules in live cells. The class of processes that we study involves the coexistence of a fractional process with correlated increments and power-law distributed sojourn immobilization times. Beyond the motion of proteins, which was studied here in detail, these processes are encountered in vastly diverse scientific fields, such as hydrology^{45,67} and movement ecology^{68}, and thus our results are expected to be widely applicable. The PSD analysis is very robust, particularly in noisy systems where it is impossible to obtain a very large number of experimental trajectories. Thus, the analysis is useful in elucidating the statistical properties of trajectories obtained by single-particle tracking in living cells, opening a new avenue in the analysis of protein transport.

## Methods

### Numerical simulations

We performed all simulations in MATLAB. To generate a CTRW (*H* = 1/2), we synthesized increments drawn from a standard normal random variable, i.e., Δ*x*^{2} = 1/2. Subsequently, the times between steps were drawn from a Pareto distribution *ψ*(*t*) = *α**t*^{−(1+α)} for *t* ≥ 1. For the subordinated random walk with *H* ≠ 1/2, we obtained the increments using the MATLAB function wfbm to generate fBM. In this case, Δ*x* is a constant that depends on *H*. For each case, a total number of 10,000 realizations were obtained with either *t*_{m} = 2^{16} or *t*_{m} = 2^{18} and a sampling time of 1.

### Live cell imaging and single-molecule tracking

Experimental details for cell culture, transfection, labeling, and imaging have been published previously^{63}. Briefly, E18 rat hippocampal neurons were plated on glass-bottom dishes that were coated with poly-l-lysine. Neurons were grown in Neurobasal medium (Gibco/Thermo Fisher Scientific, Waltham, MA, USA) with penicillin/streptomycin antibiotics (Cellgro/Mediatech, Inc., Manassas, VA, USA), GlutaMAX (Gibco/Thermo Fisher Scientific, Waltham, MA, USA), and NeuroCult SM1 neuronal supplement (STEMCELL Technologies, Vancouver, BC, Canada). For imaging, the cultures were incubated in imaging saline consisting of 126 mM NaCl, 4.7 mM KCl, 2.5 mM CaCl_{2}, 0.6 mM MgSO_{4}, 0.15 mM NaH_{2}PO_{4}, 0.1 mM ascorbic acid, 8 mM glucose, and 20 mM HEPES (pH 7.4). Neurons were transfected with a Nav1.6 construct containing an extracellular biotin acceptor domain (Nav1.6-BAD,^{63}), using Lipofectamine 2000 (Invitrogen, Life Technologies, Grand Island, NY, USA). pSec-BirA (bacterial biotin ligase) was co-transfected to biotinylate the channel. Labeling of surface channels was performed before imaging. Neurons were rinsed with imaging saline and then incubated for 10 min at 37 °C with streptravidin-conjugated CF640R (Biotium, Hayward, CA, USA) diluted 1:1000 in imaging saline. Total internal reflection fluorescence images were acquired at 20 frames/s using the 647 nm laser line of a Nikon Eclipse Ti fluorescence microscope equipped with a Perfect-Focus system, an Andor iXon EMCCD DU-897 camera, and a Plan Apo TIRF 100×, NA 1.49 objective. Imaging was performed at 37 °C using a heated stage and objective heater. Nav trajectories were obtained by single-molecule tracking using the U-track algorithm^{69}.

### Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

The datasets generated during the current study have been deposited in the Zenodo.org database under https://doi.org/10.5281/zenodo.5528301.

## References

- 1.
Mandelbrot, B. B., The Fractal Geometry of Nature (WH Freeman, 1982).

- 2.
Hu, X. et al. The dynamics of single protein molecules is non-equilibrium and self-similar over thirteen decades in time.

*Nat. Phys.***12**, 171 (2016). - 3.
Mandelbrot, B. B.

*Gaussian Self-Affinity and Fractals: Globality, the Earth, 1/f Noise, and R/S***8**(Springer Science & Business Media, 2002). - 4.
Lowen, S. B. & Teich, M. C. Fractal renewal processes generate 1/

*f*noise.*Phys. Rev. E***47**, 992 (1993). - 5.
Watkins, N. W. Mandelbrot’s 1/

*f*fractional renewal models of 1963–67: the non-ergodic missing link between change points and long range dependence. In*International Work-Conference on Time Series Analysis*197–208 (Springer, 2016). - 6.
Niemann, M., Kantz, H. & Barkai, E. Fluctuations of 1/

*f*noise and the low-frequency cutoff paradox.*Phys. Rev. Lett.***110**, 140603 (2013). - 7.
Takeuchi, K. A. 1/

*f*^{α}power spectrum in the Kardar-Parisi-Zhang universality class.*J. Phys. A***50**, 264006 (2017). - 8.
Balandin, A. A. Low-frequency 1/

*f*noise in graphene devices.*Nat. Nanotechnol.***8**, 549 (2013). - 9.
Krapf, D. Nonergodicity in nanoscale electrodes.

*Phys. Chem. Chem. Phys.***15**, 459 (2013). - 10.
Csabai, I. 1/

*f*noise in computer network traffic.*J. Phys. A***27**, L417 (1994). - 11.
Sornette, A. & Sornette, D. Self-organized criticality and earthquakes.

*EPL***9**, 197 (1989). - 12.
Ivanov, P. C. et al. From 1/

*f*noise to multifractal cascades in heartbeat dynamics.*Chaos***11**, 641 (2001). - 13.
Voss, R. F. Evolution of long-range fractal correlations and 1/

*f*noise in DNA base sequences.*Phys. Rev. Lett.***68**, 3805 (1992). - 14.
Moon, W., Agarwal, S. & Wettlaufer, J. S. Intrinsic pink-noise multidecadal global climate dynamics mode.

*Phys. Rev. Lett.***121**, 108701 (2018). - 15.
Halley, J. M. & Inchausti, P. The increasing importance of 1/

*f*-noises as models of ecological variability.*Fluct. Noise Lett.***4**, R1 (2004). - 16.
Mandelbrot, B. Some noises with 1/

*f*spectrum, a bridge between direct current and white noise.*IEEE Trans. Inf. Theory***13**, 289 (1967). - 17.
Bouchaud, J.-P., Cugliandolo, L., Kurchan, J. & Mézard, M. Spin glasses and random fields.

*Direct. Condens. Matter Phys.***12**, 443 (1998). - 18.
Sadegh, S., Barkai, E. & Krapf, D. 1/

*f*noise for intermittent quantum dots exhibits non-stationarity and critical exponents.*New J. Phys.***16**, 113054 (2014). - 19.
Rodriguez, M. A., Denis-le Coarer, F. & Valle, A. 1/

*f*noise in the intensity fluctuations of vertical-cavity surface-emitting lasers subject to parallel optical injection.*Phys. Rev. E***97**, 042105 (2018). - 20.
Leibovich, N. & Barkai, E. Aging Wiener-Khinchin theorem.

*Phys. Rev. Lett.***115**, 080602 (2015). - 21.
Dechant, A. & Lutz, E. Wiener-Khinchin theorem for nonstationary scale-invariant processes.

*Phys. Rev. Lett.***115**, 080603 (2015). - 22.
Leibovich, N., Dechant, A., Lutz, E. & Barkai, E. Aging Wiener-Khinchin theorem and critical exponents of 1/

*f*^{β}noise.*Phys. Rev. E***94**, 052130 (2016). - 23.
Kubo, R., Toda, M. & Hashitsume, N.

*Statistical Physics II: Nonequilibrium Statistical Mechanics*vol. 31 (Springer Science & Business Media, 2012). - 24.
Mandelbrot, B. B. & Van Ness, J. W. Fractional Brownian motions, fractional noises and applications.

*SIAM Rev.***10**, 422 (1968). - 25.
Montroll, E. W. & Weiss, G. H. Random walks on lattices. II.

*J. Math. Phys.***6**, 167 (1965). - 26.
Barkai, E., Garini, Y. & Metzler, R. Strange kinetics of single molecules in living cells.

*Phys. Today***65**, 29 (2012). - 27.
Goldberger, A. L. et al. Fractal dynamics in physiology: alterations with disease and aging.

*Proc. Natl Acad. Sci. USA***99**, 2466 (2002). - 28.
Yang, H. et al. Protein conformational dynamics probed by single-molecule electron transfer.

*Science***302**, 262 (2003). - 29.
Krapf, D. & Metzler, R. Strange interfacial molecular dynamics.

*Phys. Today***72**, 48 (2019). - 30.
Metzler, R., Jeon, J.-H., Cherstvy, A. G. & Barkai, E. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking.

*Phys. Chem. Chem. Phys.***16**, 24128 (2014). - 31.
Manzo, C. & Garcia-Parajo, M. F. A review of progress in single particle tracking: from methods to biophysical insights.

*Rep. Prog. Phys.***78**, 124601 (2015). - 32.
Krapf, D.

*Current Topics in Membranes*vol. 75 167–207 (Elsevier, 2015). - 33.
Sabri, A., Xu, X., Krapf, D. & Weiss, M. Elucidating the origin of heterogeneous anomalous diffusion in the cytoplasm of mammalian cells.

*Phys. Rev. Lett.***125**, 058101 (2020). - 34.
Dean, D. S., Iorio, A., Marinari, E. & Oshanin, G. Sample-to-sample fluctuations of power spectrum of a random motion in a periodic Sinai model.

*Phys. Rev. E***94**, 032131 (2016). - 35.
Krapf, D. et al. Power spectral density of a single Brownian trajectory: what one can and cannot learn from it.

*New J. Phys.***20**, 023029 (2018). - 36.
Krapf, D. et al. Spectral content of a single non-Brownian trajectory.

*Phys. Rev. X***9**, 011019 (2019). - 37.
Sposini, V., Metzler, R. & Oshanin, G. Single-trajectory spectral analysis of scaled Brownian motion.

*New J. Phys.***21**, 073043 (2019). - 38.
Sposini, V., Grebenkov, D., Metzler, R., Oshanin, G. & Seno, F., Universal spectral features of different classes of random diffusivity processes.

*New J. Phys.***22**, 063056 (2020). - 39.
Grimm, M., Jeney, S. & Franosch, T. Brownian motion in a Maxwell fluid.

*Soft Matter***7**, 2076 (2011). - 40.
Bouchaud, J.-P. Weak ergodicity breaking and aging in disordered systems.

*J. Phys. I***2**, 1705 (1992). - 41.
Vollmer, J., Rondoni, L., Tayyab, M., Giberti, C. & Mejía-Monasterio, C. Displacement autocorrelation functions for strong anomalous diffusion: a scaling form, universal behavior, and corrections to scaling.

*Phys. Rev. Res.***3**, 013067 (2021). - 42.
Margolin, G. & Barkai, E. Nonergodicity of blinking nanocrystals and other Lévy-walk processes.

*Phys. Rev. Lett.***94**, 080601 (2005). - 43.
Bel, G. & Barkai, E. Weak ergodicity breaking in the continuous-time random walk.

*Phys. Rev. Lett.***94**, 240602 (2005). - 44.
Scher, H. & Lax, M. Stochastic transport in a disordered solid. I. Theory.

*Phys. Rev. B***7**, 4491 (1973). - 45.
Berkowitz, B., Cortis, A., Dentz, M. & Scher, H. Modeling non-Fickian transport in geological formations as a continuous time random walk.

*Rev. Geophys.***44**, RG2003 (2006). - 46.
Muñoz-Gil, G. et al. Objective comparison of methods to decode anomalous diffusion. https://arXiv.org/2105.06766 (2021).

- 47.
Masoliver, J., Montero, M. & Weiss, G. H. Continuous-time random-walk model for financial distributions.

*Phys. Rev. E***67**, 021112 (2003). - 48.
Szymanski, J. & Weiss, M. Elucidating the origin of anomalous diffusion in crowded fluids.

*Phys. Rev. Lett.***103**, 038102 (2009). - 49.
Magdziarz, M., Weron, A., Burnecki, K. & Klafter, J. Fractional Brownian motion versus the continuous-time random walk: a simple test for subdiffusive dynamics.

*Phys. Rev. Lett.***103**, 180602 (2009). - 50.
Sadegh, S., Higgins, J. L., Mannion, P. C., Tamkun, M. M. & Krapf, D. Plasma membrane is compartmentalized by a self-similar cortical actin meshwork.

*Phys. Rev. X***7**, 011031 (2017). - 51.
Sokolov, I. Lévy flights from a continuous-time process.

*Phys. Rev. E***63**, 011104 (2000). - 52.
Dybiec, B. & Gudowska-Nowak, E. Subordinated diffusion and continuous time random walk asymptotics.

*Chaos***20**, 043129 (2010). - 53.
Weigel, A. V., Simon, B., Tamkun, M. M. & Krapf, D. Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking.

*Proc. Natl Acad. Sci. USA***108**, 6438 (2011). - 54.
Tabei, S. A. et al. Intracellular transport of insulin granules is a subordinated random walk.

*Proc. Natl Acad. Sci. USA***110**, 4911 (2013). - 55.
Mosqueira, A., Camino, P. A. & Barrantes, F. J. Cholesterol modulates acetylcholine receptor diffusion by tuning confinement sojourns and nanocluster stability.

*Sci. Rep.***8**, 1 (2018). - 56.
Etoc, F. et al. Non-specific interactions govern cytosolic diffusion of nanosized objects in mammalian cells.

*Nat. Mater.***17**, 740 (2018). - 57.
Levin, M., Bel, G. & Roichman, Y. Measurements and characterization of the dynamics of tracer particles in an actin network.

*J. Chem. Phys.***154**, 144901 (2021). - 58.
Roman-Ancheyta, R., de los Santos-Sánchez, O., Horvath, L. & Castro-Beltrán, H. M. Time-dependent spectra of a three-level atom in the presence of electron shelving.

*Phys. Rev. A***98**, 013820 (2018). - 59.
Bouchaud, J.-P. & Georges, A. Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications.

*Phys. Rep.***195**, 127 (1990). - 60.
Scher, H. Continuous time random walk (CTRW) put to work.

*Eur. Phys. J. B***90**, 1 (2017). - 61.
Klafter, J. & Sokolov, I. M.

*First Steps in Random Walks: From Tools to Applications*(Oxford University Press, 2011). - 62.
Shlesinger, M. F. Origins and applications of the Montroll-Weiss continuous time random walk.

*Eur. Phys. J. B***90**, 93 (2017). - 63.
Akin, E. J. et al. Single-molecule imaging of Nav1.6 on the surface of hippocampal neurons reveals somatic nanoclusters.

*Biophys. J.***111**, 1235 (2016). - 64.
Weron, A. et al. Ergodicity breaking on the neuronal surface emerges from random switching between diffusive states.

*Sci. Rep.***7**, 1 (2017). - 65.
Leibovich, N. & Barkai, E. 1/

*f*^{β}noise for scale-invariant processes: how long you wait matters.*Eur. Phys. J. B***90**, 1 (2017). - 66.
Meroz, Y., Sokolov, I. M. & Klafter, J. Subdiffusion of mixed origins: when ergodicity and nonergodicity coexist.

*Phys. Rev. E***81**, 010101 (2010). - 67.
Schumer, R., Benson, D. A., Meerschaert, M. M. & Baeumer, B. Fractal mobile/immobile solute transport.

*Water Resour. Res.***39**, 1296 (2003). - 68.
Vilk, O. et al. Ergodicity breaking and lack of a typical waiting time in area-restricted search of avian predators. https://arXiv.org/2101.11527 (2021).

- 69.
Jaqaman, K. et al. Robust single-particle tracking in live-cell time-lapse sequences.

*Nat. Methods***5**, 695 (2008).

## Acknowledgements

The Nav1.6 imaging was performed by Dr. Elizabeth Akin. D.K. thanks Dr. Mike Tamkun for his help with the experiments and useful discussions. We acknowledge the support from the Colorado State University Libraries Open Access Research and Scholarship Fund (D.K.), National Science Foundation grant 2102832 (D.K.), and Israel Science Foundation grant 1898/17 (E.B.).

## Author information

### Affiliations

### Contributions

D.K. conceived and designed the project, and constructed the theory. D.K and Z.R.F. generated numerical simulations and analyzed the data. D.K, E.B., and Z.R.F. interpreted the results. D.K. wrote the manuscript, assisted by E.B. and Z.R.F.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Fox, Z.R., Barkai, E. & Krapf, D. Aging power spectrum of membrane protein transport and other subordinated random walks.
*Nat Commun* **12, **6162 (2021). https://doi.org/10.1038/s41467-021-26465-8

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.