Fractal properties of ion channels and diffusion

Fractal Properties ROBERTO

RONCAGLIA,

of Ion Channels RICCARDO

and Diffusion

MANNELLA,

and PAOLO GRIGOLINI*

Dipatiimento di Fisica e GNSM de1 CNR, 56126 Pisa, Italy Received 26 May 1993; revised 6 December

1993

ABSTRACT We focus our attention on a fractal model recently proposed by Liebovitch to account for the lack of a time scale in ion channel kinetics. We establish a connection between the dwell-time distributions and the correlation time of the ion channel signal, thereby making it possible to derive analytical predictions on the diffusion properties of a random walk constructed from the sum of the current fluctuations of ion channels. With the help of a numerical simulation of the Lievovitch model, it is shown that the Hurst analysis can provide a reliable determination of the standard (or anomalous) diffusion properties. On the basis of results of computer simulation we argue that by applying the Hurst analysis to the experimental distribution of closed times it is possible, in principle, to establish whether the Liebovitch model is valid.

1.

INTRODUCTION

In the recent literature on ion channel kinetics there has been a controversy [6,7] as to whether the fractal model proposed by Liebovitch [lO,ll] or the many-state Markov model proposed in the standard analysis of the experimental recordings [2] affords a proper representation of ion channel kinetics. Liebovitch [7] pointed out that the Markov model need not be preferred to the fractal model even if the Markov model were a better fit to data than the fractal model. This is so because the better accuracy would only reflect the possibility of changing a larger number of fitting parameters.

*Istituto di Biofisica de1 CNR, 56127 Pisa, Italy, and Department University of North Texas, Denton, TX 76209.

MATHEMATICAL

BIOSCIENCES

of Physics,

123:77-lOl(1994)

OElsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010

77 0025-5564/94/$7.00

7x

ROBERTO

RONCAGLIA

ET AL.

More recently, Liebovitch [S] showed that all the statistical properties of a fractal model with only two states but time-dependent transition coefficients can be reproduced by suitable Markov models. He singled out two such examples that are able to reproduce the statistics of the bistable fractal model, each of the two resting on the assumption that the open state is coupled with an infinite number of closed states. Here the transition rate was chosen to satisfy a decreasing power law, so as to mimic, also in Markov systems, the absence of a time scale characteristic of fractal processes. In this paper we focus our attention on the bistable fractal model proposed by Liebovitch, slightly modifying the time-dependent transition probability originally introduced by him so as to give the open state a finite rather than a vanishing statistical weight. Then we evaluate analytically the correlation time of the time-dependent signal, and we show that it is characterized by a power law. This allows us to make a well-defined prediction on the anomalous or standard nature of the diffusion process associated to the dynamics of the ion channel. This is the diffusion process corresponding to a random walker that takes a step ahead by a fixed quantity for any transition from the closed state to the open state. We establish a precise connection between the power-law decay of the correlation function of the ion signal and the anomalous character of the corresponding diffusion process. Then we illustrate an efficient procedure for analyzing the experimental data to assess whether or not anomalous diffusion exists. Its existence would lend strong support to the Liebovitch model. The plan of the paper is as follows. Section 2 is devoted to illustrating a generalization of the model introduced by Liebovitch, for the main purpose of giving a finite rather than a vanishing statistical weight to the open state. Section 3 shows how to evaluate analytically the correlation function of the ion signal on the basis of the generalized Liebovitch model. The main result is that the fractal nature of the Liebovitch model results in a correlation function with power-law decay. In Section 4 we establish a connection between such a power-law decay process and the anomalous nature of the diffusion process. In Section 5 we illustrate an efficient numerical analysis to assess, on the basis of real experimental data, the expected anomalous nature of the diffusion process. In Section 6 we discuss the results of a computer simulation on a stochastic system equivalent to the generalized model of Liebovitch under study. It is shown that the method proposed for analyzing the experimental data might reveal the power-law decay of the signal correlation function even in a physical condition of normal diffusion. Section 7 is devoted to pointing out the original results of this paper and the benefits that might be derived from an analysis of the experimental data carried out according to our prescriptions.

FRACTAL

2.

79

PROPERTIES

THE FRACTAL MODEL UNDER STUDY

Experimental observation of the dynamics of ion channels seems to suggest that no well-defined time scale exists [7]. To account for this important property, Liebovitch [lo, 111 proposed a two-state model with one open state and one closed state and with time-dependent probabilities for moves between states. The model is schematically illustrated in Figure 1. More specifically, Liebovitch [lo, 111 assumed that the time dependence is exhibited by the probability that the system will move from the closed state to the open state and that this transition probability becomes smaller and smaller as the time spent in the closed state increases. The dependence of this transition probability is assumed to have a power-law form, which is then shown [lo, 111 to be responsible for the lack of a time scale in the probability distribution of the dwell times. We adopt a slightly modified version of the Liebovitch model, defined as follows. The transition probabilities k,,(t) and k,(t) are given the form

k,,(t)

= A(C + tp,

l
k,,(t)

= h.

(2.1)

Liebovitch’s original version is recovered by assuming that t, = 0. However, this would make the transition probability k, diverge at t = 0, thereby making the statistical weight of the open state vanish. The price we have to pay to make finite the statistical weight of the open state is to renounce having a complete lack of time scale. In the short-time region with t = t,, the properties of self-similarity are not yet exhibited; they appear only for t x=- t,. On the other hand, this is a quite plausible property, and it is confirmed by the numerical analysis of many dynamical systems, where the lack of a time scale is proven to make itself manifest only in the long-time region [13].

(-J-g+-) FIG. 1. Scheme of the Liebovitch model for the ion channel activity. 0 stands for open state and C for the closed state.

80

ROBERTO

RONCAGLIA

ET AL.

It is straightforward to show that the generalized version of the Liebovitch model (2.1) leads to the following open and closed dwell-time distributions:

fdt) = hexp(- ht) A[t + t,,y

(2.2a) exp[&~~~“]exp[-&(~+,,,)2~d],

i

f,(t) =

l
+ 4,1’+ n, d=2.

(2.2b)

It must be pointed out that only the case d = 2 results in a closed time distribution with power-law decay property. It is also worthwhile to notice that, according to Liebovitch [lO,ll], precisely this value of the fractal exponent d satisfactorily fits the corresponding experimental data. The present paper is essentially devoted to the analysis of this interesting case. Within this context, our major purpose is to relate the property of distributions (2.2) to another statistical property, namely, diffusion. We shall prove that the fractal behavior in the distribution of dwell times might result in anomalous diffusion. To prove this interesting connection between the fractal nature of the long-time distribution of dwell times and anomalous diffusion, we must evaluate the correlation function of a dichotomous variable 5, which defines the state of the ion channel, and its n-time joint or conditional probabilities. The theoretical predictions on the dwell-time distributions are usually the final goal of a theoretical analysis. We follow a different avenue. We rely on the main conclusions of Liebovitch [7,10] and assume that the statistical properties of Equations (2.2) fit satisfactorily the corresponding experimental data. Thus Equations (2.2) are the starting point of our theoretical investigation. We also assume that the statistical properties of the ion signal are completely defined by the dwell-time distributions (2.2). This assumption surely holds true in the model we are considering. In fact, it is easy to realize that the probability of having a closed (open) interval of length t is totally independent of the previous sequence of intervals and that, in other words, the dichotomous signal generated by the Liebovitch model is statistically indistinguishable from a signal obtained alternating open and closed intervals distributed according to distributions (2.2). It must be stressed, however, that although in the model here under discussion this property holds true, and it turns out to be beneficial, it is not necessarily satisfied by real systems [Il.

Xl

FRACTAL PROPERTIES

3.

CORRELATION

The purpose of the statistical

FUNCTION

of this section is to determine the correlation variable 5, with the two values

tc=o,

t,=c,

function

(3.1)

where 5, denotes the value of the variable 5 corresponding to closed state and .$(, refers to the open state. The value of this statistical variable corresponding to the closed state is set equal to zero for the sake of calculational convenience. The correlation function is given by definition by the expression

where P,(t,,ti) is the probability for the stochastic variable 5 being t, at time t,, while P2( &, t,l cl, t,> denotes the conditional probability, that is, the probability that 5 will have a value .$ at time t, given that at time t, it was tj. Since in this paper we consider only stationary processes, we can replace P,( tj, tl) in Equation (3.2a) with the time-independent equilibrium probability, which we denote as P,( tj,m). In conclusion, using prescriptions (3.0, we can rewrite the expression for the correlation function only in terms of the open-open conditional probability,

C(f) = wbw)) = cP,(open,t

3.1.

ANALYTICAL

EVALUATION

= ~2q~“,~Itooo)P~( 6”,M) lopen,O).

OF THE CONDITIONAL

(3.2b)

PROBABILITIES

Before we evaluate C(t) we must determine the conditional probability. The stochastic process under study is an alternating renewal process. It is straightforward to realize that, since the transition probabilities (2.1) depend on the time spent in the closed state, it is not possible to write down directly a master equation for the conditional probability. The calculation must therefore follow a different strategy. Let us start looking for an expression for the conditional probability and in particular for P, (open,rlopen,O)-that is, the probability that the system, originally in the open state, will again be in the same configuration at a given time t. This can be calculated using the following path-integral-like technique. Let us consider all the possible sequences

ROBERTO

=,

RONCAGLIA

ET AL.

o
t

I

I

t1

t2

t

I I t1

t3

tr

t

I

I,

t-l

t

0 -c t1-c t2<

0 <

t

t1< tz < t3< t‘4-ct

I

L t1

t2

t3

--I

0 < t1< t2-c . . .. -z t2n-I < t2n-c t2n-2

t2n-I

C2n

t

t

FIG. 2. All possible paths that can be followed from the open level at t = 0 back to the open state at time t. The upper level corresponds to the open configuration of the channel, the lower to the closed one.

of open-closed intervals that, starting from an open configuration after a time t, end up in another open state. Following the scheme of Figure 2, we can label the different paths satisfying the above condition by counting the number of closed states that occur in the considered interval of time. Thus set 0 contains all the paths that never flip down to the closed state, set 1 contains those that flip only once, and so on. Let us now denote by P(n, t) the probability for set it, that is, the sum of the probabilities of the individual paths contained in set II. Since all the paths are statistically independent of each other, we can write the conditional probability as

P,(open,tlopen,O)

= f P(n,t). n=f.l

(3.3)

We now have to determine the probabilities P(n, t). The first one, P(0, t), is easily found by noticing that the paths contained in set 0 are those with open intervals longer that t, that is,

(3.4)

83

FRACTAL PROPERTIES

where f,(t) is the open-time distribution. The calculation of the next probability P(l,t), implies the sum over all possible lengths of the closed interval and over all possible lengths of the final open intervals. This gives the expression

where f,(t) is the closed-time distribution. Applying the same arguments, it is straightforward to generalize the above result and obtain for P(n, t) the expression

By a suitable change of integration variables, Equation transformed into the product of convolutions

P(n,t) =

(3.6) can be

[f,(t).fc(t).f,(t).f,(t)...S,(t).f,(t)l.~(t), (3.7)

where the term in brackets is equal to equivalent to

P(n,t) =f,(f)-f,(t)-P(~

[fJt>*f$)l”. Equation (3.7) is

-l,t),

n = 1,2,3 ,....

(3.8)

Note that in Equations (3.7) and (3.8) we have used the definition

f(t).g(t) =

pe - T)dT).

(3.9)

Equation (3.8) serves the useful purpose of providing a compact expression for the conditional probability. Indeed, by inserting (3.8) into (3.3), we obtain the following integral equation for P, (open,t lopenO): P2(open,tlopen,O)

=~(t)+f,(t)~fC(t).Pz(open,tIopen,O).

(3.10)

84

ROBERTO

RONCAGLIA

ET AL.

Using the well-known property of the Laplace transform of a convolution product, the Laplace transform of the solution of (3.10) is easily written as

P,(open,zlopen,O)

.6(z)

= l-

(3.11)

f,,(Z)SC(Z)

.

The other conditional probabilities, P,(closed, t 1closed, O), etc., are easily found following the same arguments as those used earlier. In particular, P,(closed, t (closed,01 is obtained from Equation (3.11) by switching the open distribution with the closed one, thereby obtaining

P,(closed,z

Iclosed,O)

L(z)

=

l-f,,(Z)fC(Z) The cross-conditional terms of convolution

(3.12) .

probabilities are obtained by expressing them in products of the previous probabilities as follows:

P,(open,rlclosed,O)

=/dTfC(7)Pz(open,t t

-T/open,O),

P,(closed,tjopen,O)

=ij’drfO(r)Pz(closed.t

-Tlclosed,O).

(3.13a) (3.13b)

All the n-point joint or conditional probability functions can be obtained by applying the same technique. However, since we are interested only in the statistical properties of the model under study at one and two times, we shall not derive their expressions. In the remaining part of this paper it will turn out to be useful to study the stationary properties of the occupation probability P,(t, t>. By definition this is readily obtained by replacing Equations (3.11)-(3.13) in the expressions P,(closed,r)

= P,(closed,r + P,(closed,t

P,(open,r)

CORRELATION TIES

lopen,O)P,(open,O),

(3.14)

= P,(open,tlopen,O)P,(open,O) + Pz(open,t

B.

Iclosed,O)P,(closed,O)

FUNCTION

Iclosed,O)P,(closed,O).

WITH d = 2 AND

Liebovitch has shown [9] that in the channels the main statistical characteristics

ITS ASYMPTOTIC

PROPER-

case of K + selective ion of the signal, such as the

FRACTAL PROPERTIES

85

dwell-time distributions, are well fitted by the fractal model (2.1) provided that the fractal dimension d is equal to 2. It is worthwhile to note that in Liebovitch’s original work the transition probability K,,(t) is simply set to be t’ -d so that the experimental value d = 2 would imply an unnormalizable closed-time distribution. This singularity is essentially due to the divergent character of K,,(t) at t = 0. Actually, it is clear that, rather than being a genuine property of the experimental recording, this feature derives from the particular form of K,,(t) chosen to describe the system. Nevertheless, the experimental data show that for sufficiently long time, the transition probability does scale according to the power-law function proposed by Liebovitch. We are thus led to the conclusion that the fractal property of the dwell-time distribution is a time asymptotic property of the system. By the same token we are led to expect that all the manifestations of the fractal nature of the process will take place only in the asymptotic regime. In the model we are considering, this corresponds to times t z+ t,, [see Eq. (2.111. Let us therefore try to gain information on the fractal nature of the process by studying the asymptotic regime of the correlation function of our system. Before doing this it is important to exclude from the analysis those values of the parameters in Equation (2.1) that lead to nonphysical situations. To this purpose let us calculate the stationary probabilities of our system, that is, the limit t -00 of Equation (3.14). This result is easily obtained by evaluating the limit z + O+ of the Laplace transform of (3.14) and reads

0, P,(open,m)

=

for A G 1,

A-l A - 1 + ht, ’

forA>l.

(3.15)

We consider that the condition with a vanishing statistical weight for the open state is not quite correct on physical grounds. Therefore, we do not consider the case A G 1. Let us now focus our attention on the time evolution of P,(open, t (open,O), which, as we have seen, is proportional to the correlation function of the signal. Using in Equation (3.11) the explicit form of fO(z>, the Laplace transform of P,(open, t lopen,O) becomes

P,(open,

z Jopen,O)

= 1 z ( l+h&z)).

(3.16)

86

ROBERTO

RONCAGLIA

ET AL.

where P,(z) is the Laplace transform of the probability of Z’,(t) remaining in the closed state during the interval of time [O,t] without making any transition to the open state,

In the case we are considering

[d = 2 in Eq. (2.2)1, this reads

fc( z) = t,;le-z’,lZ A-‘r(lwhere

IX1 - A,t,z)

is the incomplete

A,t,,z),

gamma

function

(3.18) [4] defined

cc ( -l)n(t”Z)‘-A+n c rt!(l-A+n) n=O

I(l-/I,t,z)=T(l-A)+

A#1,2,3

as

’ (3.19)

,.a..

Note that for the sake of simplicity we omit the discussion of the case of integer values of A, which, although leading to a different expression of I on (3.181, does not alter the general properties of the asymptotic behavior of P,(open, t Iopen, 0). To evaluate the asymptotic behavior of the correlation function we shall make use of the Tauberian theorems, and hence we need to determine the behavior of the Laplace transform (3.16) around z = 0’. To this end it is useful to rewrite P,(z) as (3.20)

Pc(z)=t,/(A-l)+~c(z). From Equation

(3.18) it is easily found ,ty+

thereby around

E(z)qT(l-

that A)zA-‘,

(3.21)

suggesting a straightforward way of expanding Equation z = O+. In fact, using Equation (3.20) this becomes

Pz(open,

z /open,O)

(3.16)

=

p,“(z) - ...

1

, (3.22)

FRACTAL PROPERTIES

87

which, together with (3.20, gives as leading asymptotic terms for z = O+ the expression P2(open, 2 lopen,O) = i

1

0 i fle+‘,, A

!

-h

t,“r(l-

A)z~-~.

(3.23)

Finally, using Equations (3.2b) and (3.151, the asymptotic behavior of the correlation function is obtained by means of a Tauberian theorem and reads

l+&:I:ht(1 2

A-l

0

1

,

t

z+ to.

(3.24)

The power law in Equation (3.24) is a clear fingerprint of the fractal nature of the statistical process generated by the bistable model (2.1). The correlation C(t) is a self-similar function of time, thereby implying a decaying behavior with no intrinsic time scale (obviously, if we neglect the transient regime, t = to>. In the next section we shall see that this power-law behavior can affect the diffusive properties of the stochastic signal, resulting in an anomalous diffusion, and that an efficient way of detecting anomalous diffusion from “experimental” data is provided by the Hurst analysis. We shall prove this important aspect by applying this procedure to the results of numerical simulation rather than to real experimental data. Before ending this section, it is convenient to explain the meaning of diffusion within the context of the dynamics of ion channels. By diffusion properties we mean the asymptotic behavior of the second moment of the statistical time-dependent process x(t) constructed from the zero-average statistical variable i(t) = t(t) - ( .$)eq by adding up all the values of g in the interval of time 10,t], x(t)

‘pTg(7).

(3.25)

In the case of vanishing correlation time, when the correlation function reaches its equilibrium value instantaneously, with the mere help of the central limit theorem it is easy to show that the stochastic variable x(t) has a “normal” diffusive behavior, that is, the second moment of x(t) increases linearly with time.

(x”(t)> = Dt.

(3.26)

ROBERTO

88

RONCAGLIA

ET AL.

This is the conventional diffusion process, and it implies that the correlation function C(t) relaxes very quickly, if not instantaneously, to its equilibrium value. Actually, the generalized Liebovitch model that is adopted in this paper makes the correlation function C(t) result in the time-asymptotic property of Equation (3.24). This is not a quickly relaxing correlation function, and the consequences of this slow decay of C(t) on diffusion will be studied in the next section. 4.

DIFFUSION Let us consider

REGIME the dynamical

equation

where z(t) denotes the time-dependent statistical process 2 = t(t)( 5 )eq. We have in mind the dichotomous stochastic variable c(t) of Section 3. Ths does not have a vanishing mean value, and for this reason its correlation function does not relax to zero for t tending to x. The variable 2 on the contrary, fulfills the property

( 8 Lq = 0, and its normalized

correlation

function,

a(t),

(4.2a) defined

by

@(t) = (5(0)5(w?_q/(me,, is related

to the correlation

Q(t)

function

=

(4.2b)

C(t) by the expression

C(t)-w$ ( 5 2)eq- ( 5 2, .

(4.3)

The subscript eq means that the averages are carried out on the stochastic process 5, which is supposed to be stationary. It seems to be quite plausible to apply this assumption to the ion channel signal. Notice that Equation (3.25) is the solution of (4.1) corresponding to the initial condition x(0) = 0. By evaluating the square of (3.25) and taking the average over many random trajectories corresponding to the same initial condition, we get

(4.4)

89

FRACTAL PROPERTIES

By making stationary,

the reasonable namely that (

Equation

assumption

,p)~(t”)>c,

=

that the stochastic

(5(O)@‘-

process

w-v

g(t) is

(4.5)

(4.4) is shown to lead to

(4.6) The standard diffusion of (3.26) is recovered if the correlation function Q,(t) decays to zero quickly enough to make finite the correlation time 7=

/0

%(l)

dt.

(4.7)

In this case for times t of the order of T z- 7, the upper time integration limit on the right-hand side of (4.6) can be replaced by infinity, and the integration of the resulting equation leads to (3.26) with D=($2),q~. Let us consider now the case where the correlation function is given by Q(t)

---)

(4.8) the time-asymptotic

behavior

(4.9)

Ktpp.

From (3.24) we see that the ion channel model we are using precisely to a correlation function of this kind, with PEA-l. We apply the further

which implies

condition

of

leads

(4.10)

that

P >o,

(4.11)

A>1

(4.12)

the condition

and consequently [see Eq. (3.1511.

a nonvanishing

statistical

weight

for the open

state

90

ROBERTO

RONCAGLIA

ET AL.

In this new condition, it is expected that the standard diffusion (3.26) is no longer valid and that the long-time regime is described

(x2(t)) = Dti,

of by

(4.13)

with 5 not necessarily equal to 1. In the specific case where p is large enough as to make the correlation time r finite, in the long-time limit we must again recover the standard diffusion process. Thus we have l=l In the cases i > 1, differentiating (4.13) result must coincide once with respect to

if /3>1.

(4.14)

we establish a connection between < and p by twice with respect to time and noticing that the with the equation obtained by differentiating (4.6) time. We thus have l=2--/3

if 1 2 p > 0.

The reverse is also true. If we have anomalous immediately conclude that the correlation function power law and that p=2-5.

(4.15) diffusion, then we Q,(t) decays with a

(4.16)

Of course, in the case of standard diffusion, in principle we cannot distinguish the case with /3 > 1 from that where the correlation function undergoes another kind of fast decay, for instance, exponential decay. For this reason the experimental detection of normal diffusion would make it impossible to reach an easy conclusion as to the relevance of the Liebovitch model. However, in Section 6 we propose a technique of analyzing experimental data (by applying it to the results of numerical simulation) that makes it possible to reveal power-law decay, if it occurs, even in the case where the asymptotic time limit is characterized by normal diffusion. This technique, virtually equivalent, in the case of the model here under study, to the determination of the second moment (x2(t)> but more efficient than this method, allows us to dig out from the process of standard diffusion the existence of a power-law behavior for the correlation function Q(t). We have seen that Equation (4.6) is an exact equation that does not need any special condition. However, (4.6), through (4.14) and (4.15), only provides information on the time scale properties of ( .x”(t))“2. We do not know yet if (.x~‘Y~))“(~~) scales with the same law. Using the projection operator method and a perturbation analysis rigorously kept at the second order [5], it is possible to find the equation of motion

91

FRACTAL PROPERTIES

driving the time evolution of the probability distribution of the variable p(x,t). This equation reads (4.17) It is straightforward to prove that this equation leads to the same time scaling for all moments. This is so because (4.17) fulfills the Gaussian property (p(t))=(2n-1)!!(ap(t))n.

The diffusion equation (4.17) is diffusion coefficient. However, as totic time limit Equation (4.17) is equation by a proper time scaling, by t* = tL. 5.

(4.18)

characterized by a time-dependent already noted in [13], in the asympturned into a conventional diffusion namely, by using the time t* defined

HURST ANALYSIS

Let us consider a time-dependent average over the interval of time r: W,=+

signal

and define its time

k 6(t).

(5.1)

f=l

Let X(t,T) be the accumulated departure of t(t) from the mean ( 07,

X(t,7) = i

{t++W7}.

(5.2)

u=l

Introduce

M(T), defined by

M(T) = lr_yxTX(t,r)

. .

- lm;lqTX(t,7), d

.

(5.3)

where t is some discrete integer-valued time and T is the time span considered. Let us consider also the standard expressin for the variance,

w=; i [s(t)-m]* i t=1

10

I.

(5.4)

The Hurst analysis 131consists in studying the ratio R(T) = M(T)/V(T).

(5.5)

92

ROBERTO RONCAGLIA ET AL.

Hurst found that for many records empirical relation R(r)

in time, R(T) is well described

arH.

by the

(5.6)

In the case when the fluctuation t(t) is that of the fractional Brownian motion, Mandelbrot 131 (see also 1121) proved that the exponent i of (4.13) is related to the Hurst exponent H by the simple relation H = c/2.

(5.7)

By time scaling arguments it is shown [13] that (4.17) also in the limit of large times must become equivalent to the fractional Brownian motion. Thus, on the basis of Mandelbrot theory we conclude that any stochastic or chaotic process that can be described by the diffusion equation (4.17) in the long-time limit should be characterized by a well-defined Hurst exponent. Note that according to our analysis, the correlation between the asymptotic Hurst exponent and the power p of the correlation function a(t) should be H = 0.5

ifp>l

(5.8)

if 1 2 /? > 0.

(5.9)

and H=l-p/2

Of course, under the assumption that the correlation function Q,(t) in the long-time limit follows a power law, it is possible to derive p from H. If H = 0.5, p must be larger than 1, and for 0.5 < H < 1 we get p=2-2H.

(5.10)

Note that the case H = 0.5 does not necessarily imply a correlation function Q(t) with a power-law decay. Thus, a nonfractal channel-for example, one with single exponentially distributed open and closed times-would result in H = 0.5. It must be noted that if Equation (4.17) is true, then the Hurst coefficient can be identified with the exponent 5 appearing in (4.13). In fact, with a proper time scaling, that is, using t” = ti, Equation (4.17) is changed into a new one that in the long-time limit becomes the standard diffusion equation, leading to a Hurst coefficient H = 0.5. Coming back to the original time, this yields H = l/2 [13]. Using more direct arguments, we can say that since the Hurst coefficient refers to and in the Gaussian case all the moments the time scaling of (P)“*,

93

FRACI’AL PROPERTIES

scale in the same way, then coefficient of (x2)“*.

H

6.

OF THE FRACTAL

NUMERICAL

ANALYSIS

can

be identified

with

the

scaling

MODEL

To dig out a power-law behavior from the region at the border with normal diffusion we must carry out an analysis of the time evolution of the second moment more accurate than that which led us to (4.14) and (4.15). We assume that the long-time behavior of the correlation function @(t) is given by (4.9). However, this correlation function must also fulfill the initial condition Q(O) = 1. This means that in the short-time region the correlation function must depart significantly from the power-law behavior, which would give a divergence for t + cc).In other words, for times on the order of a certain time rc, the decay of the correlation function undergoes a decay process much different from the power-law behavior, for instance, an exponential decay process. For times larger than rc, the correlation function decays with a power law. Let us discuss the consequences of the deviation in the short-time region of the correlation function Q(t) from a power law. Notice that Equation (4.4) can also be written as (6.1) Consider the large time limit of the correlation function Q(t), that is, the temporal region with t zs=-TV. Dividing the range of integration into two subintervals [O,T,] and [r,, t’], the result of the first integration in (6.1) yields

Ql( t’) = i


dt’),

t’

K, + K,PP,

t’ >

Tc,

(6.4

where g(t’) comes from the integration of the correlation function in the time region where the asymptotic power law does not hold. Substitute (6.2) into (6.1) and consider the time region t s==TV.The integration over the first interval of time, namely [O,T,], contributes to the final results only with a constant factor. The leading contributions to (x’(t)> derive from the second interval of integration: the constant term on the right-hand side of (6.2) results in a linear increase of the second moment, whereas the second term agrees with the analysis of Section 4 and results in a power-law increase of the second moment, with exponent 2p. In conclusion, the long-time limit of (6.1) has the general form

=

[c,+C,t

+ c,t*-P])

t >> 7c.

(6.3)

94

ROBERTO

RONCAGLIA

ET AL.

This means that in the asymptotically large time region the third or second term on the right-hand side of (6.3) dominates, according to whether the condition 1 < A < 2 or 2 < A < 3 applies [for A > 3 the last term of (6.3) vanishes for t -aI. Using Equations (5.8) and (5.9) we are led to the following theoretical prediction: H=0.5

forA>

(p>l)

(6.4a)

and H=lS-

A/2

forl
(l>p>O).

(6.4b)

Thus, we see from Equation (6.3) that upon an increase in time, either ordinary or anomalous diffusion tends to become predominant according to whether our theoretical analysis of Section 4 results in the standard (H = 0.5) or in the anomalous (H > 0.5) Hurst coefficient. However, around the condition A = 2, the ordinary and anomalous contributions to (6.3) are of comparable importance (unless the extremely long time region is examined). Thus the analysis of the results of numerical simulation requires a method of high accuracy. Note that the Hurst analysis provides the same information as that already contained in the correlation function; the main advantage is that it minimizes the numerical error associated with the statistics for the same amount of CPU time. We can give some intuitive explanation for this. First of all, note that typical correlation function will decay (perhaps with some oscillations) as time increases; this implies that the power-law tail will, in many cases, set in at times for which the correlation function itself is very small, but in analyzing real data the limited statistics will very often make impossible a direct analysis of a very small correlation function. Also, it is often very important to be able to work out properties related to the correlation function over a long time domain. In this case the direct calculation of the correlation function becomes expensive, in terms of CPU time. If N is the number of data items to be analyzed, a fast Fourier transform technique (the fastest algorithm to compute correlation functions) will require a number of operations that will increase as N In N, whereas the Hurst analysis will require only a number of operations that goes as N. One could think of using the spectral density at short frequencies to derive the behavior at large times; again, this requires a number of operations that increases (at best) as N In N. Note that we are assuming that, for the correlation function or the spectral density, the (time or frequency) intervals at which the quantity is computed are evenly distributed so that we are able to apply an FFT technique. Obviously, the typical

FRACTAL PROPERTIES

95

power-law dependency of the anomalous diffusion is best studied on a log-log chart; this implies that most numerics done using either the correlation function or the spectral density, referring to the short time domain, wiil be useless. Trying to compute the correlation function or the spectral density at logarithmically spaced intervals makes the length of the calculation proportional to N’. Similarly, the direct calculation of the time-dependent moments requires N 2 operations. To make more transparent the difficulties associated with the analysis of real experimental data, we carried out a numerical experiment. The computer experiment serves a twofold purpose. First it leads us to derive an efficient method of analysis for the results of computer simulation, which would turn out to be especially important when dealing with a real experiment. Second, it has the purpose of validating the reliability of the Hurst analysis of Section 5 to assess the anomalous character of the diffusion process. Using a stochastic algorithm we simulated the two-state model of Section 3 with an equivalent stochastic process. The distribution of residence times in the open state was obtained from straight Monte Carlo generation of an exponential function. For the distribution of residence times in the closed state, we tried two procedures; we either wrote a stochastic differential equation for the escape from the closed state, which, integrated, yielded the distribution of exit times, or, alternatively, we theoretically worked out the distribution and generated the exit times accordingly. The former method suffers from an intrinsic slowness; to generate the exit times with reasonable accuracy, the integration time step of the stochastic differential equation must be fairly small, equal to t, /lo, which leads to the generation of a large number of random deviates, with excessive CPU time expense. The latter method, though, has some drawbacks also; it is interesting to appreciate that, particularly when A is close to 1, the distribution of exit times is characterized by a fairly long tail (see Figure 3); for the A’s considered (see Figure 4) it is necessary to generate exit times up to 10’ times the closed-state lifetime. A straight rejection method would be far from efficient given the small portion of area under the curve compared to the total area when large times are considered. Thus, we eventually used a mixed strategy, which consists in exactly mapping the exit time distribution up some large time (on the order of lo3 times the residence time in the closed state) onto a different probability distribution for which a rejection method is more efficient. If, on the contrary, the exit time had to be generated in the tail, then we used the stochastic differential equation. In Figure 3 we compare the dwell-time distribution resulting from the numerical experiment and the theoretical prediction, and the agreement looks very satisfactory.

96

ROBERTO

RONCAGLIA

ET AL.

1 0-l 1 o-2 1 o.3

C

va 1 o.4 1 o.5 1 o-6 1 oe7 0

100

200

300 t

400

500

FIG. 3. Probability distribution of the dwell times of the closed state with A = 1.3. The squares indicate the results of numerical simulation, and the full line indicates theoretical prediction.

The sequence of residence times in the open and closed states was then expressed in terms of discrete times, arbitrarily assigning (see Section 3) the value 0 to the closed state and 1 to the closed state, up a large time T,,,. We typically used 5 x 10’ bins and a binning time of approximately one-tenth of the residence in the closed state. We had always t, = 0.1. We then applied the analysis of Section 5 to this sequence. To improve the statistics of the numerical experiment, in the case of T

0.1 1

10

100

1000

‘s FIG. 4. The Hurst coefficient

R(T) for various values of A.

10000

FFUCTAL

PROPERTIES

97

(much) smaller than T,,,,, we divided the random sequence in the subset of sequences of length r, and we averaged the R’s [see Eq. (5.511 obtained by applying the Hurst analysis to each subset. Note also that given the power-law dependence of R(T) on 7, we took the T’S distributed exponentially; for all but the largest r’s (see Figure 4) the statistics of our computed R(T) is extremely good. Finally, 100 values of R(T) were averaged for each value of A to obtain the curves shown in Figures 4 and 5. The results of the simulations are summarized in Figures 4 and 5. Note that the R(T) curves are extremely smooth. Figure 4 shows the quantity R(T) vs T: the power-law dependency is very clear. For a standard diffusion process we should have that R(T) - TV/'. A simple way to check whether our system undergoes standard diffusion is to look at the behavior of R(T)/T 1/2 for large T. For standard diffusion we should have a horizontal straight line. The quantity R(T)/T~/~ is plotted in Figure 5; we see that the diffusion is likely to be nonstandard. As shown previously, however, we do not expect standard diffusion for the fractal model but rather a more complex functional dependency of R(T) as function of T [Eq. (6.3)]. To check in detail whether Equation (6.3) correctly describes the diffusion process, we decided to fit the quantity R(T) obtained from the numerical simulations. Again, it is not straightforward to carry out the fit; it is necessary, in fact, to decide which weight is the most appropriate given the exponential variation of the data value over the range of T considered. In general we decided to weight the data with l/r. Although this may look rather peculiar, we

r\-

I

.”

I I .. . .

. . ..

A=,

.J

----- A~l.7 --A=1.9

FIG. 5. The Hurst coefficient

a horizontal

straight line.

divided by 7 I/‘. Standard diffusion should result in

98

ROBERTO

RONCAGLIA

ET AL.

point out that even had we chosen slightly different power law for T we would have obtained consistent results for H. There is no doubt that we need some kind of decreasing weighting function; otherwise we would end up fitting properly only the data at the largest T value we considered. The fitting function also must be chosen appropriately. Table 1 illustrates the results of two different fitting procedures. The first was carried out from the mere analysis of the long-time asymptotic limit, namely, fitting the data of computer simulation with R(r) proportional to TV. The results are illustrated in the second column of Table 1. The second fitting procedure was carried out as follows. First the square root of Equation (6.3) was identified with R(t). Then the parameters C,, C,, C,, and p were used as free fitting parameters to get the best agreement with the data of the computer simulation of Figure 4. The value of H was obtained as the asymptotic limit of ln[ R(r)]/(ln 7). All the values of the third column of Table 1 but the last one were obtained using this procedure. We see that the result of the analysis of the data of computer simulation leads to values that are systematically closer to the theoretical prediction than the values resulting from the former fitting procedure.

TABLE I H

A

First Fit”

Second Fit h

1.3 1.5 1.7 1.9 2.1

0.803 0.712 0.638 0.594 0.558

0.805 0.721 0.650 0.573 0.472*

Theoretical

Value

0.850 0.750 0.650 0.550 0.500

“The Hurst exponent is obtained by fitting the “experimental” data of Figure 4 with the analytical function R(T) = 7”. hThe experimental data of Figure 4 are fitted using the square root of Equation (6.3). The last value in this column (denoted by* I is the value l- p/2, with /I obtained from the best fitting between the square root of Equation (6.3) and the “experimental” data of Figure 4. It suggests the possibility of detecting a power-law decay of the correlation function C(t) even when the Hurst coefficient H = 0.5 would not make it possible to distinguish a power-law decay from the exponential decay of the correlation function. ‘Theoretical prediction. H is given by the formula H = (3 ~ A)/2 for A < 2 and H = 0.5 for A > 2.

FRACI-AL PROPERTIES

99

The value reported in the last row of the third column is of special interest. This was obtained from the second fitting procedure, identifying H with 1- p/2, with the value of p provided by the fitting procedure itself. The theoretical value of j3 would in this case be 1.1, and the analysis of the data of computer simulation would lead in this case to p = 1.056, thereby suggesting the possibility of revealing the presence of a power-law decay even in the region where the Hurst coefficient H = 0.5 would not make it possible to distinguish a quickly decaying correlation function (for instance, a conventional exponential decay) from a power-law decay. We make a few remarks about the overall agreement between computed and theoretical H values. The table suggests that this agreement is no worse than some 5-7% (see the entry relative to A = 1.3, the one with the largest discrepancy). We regard this as a very satisfactory agreement; the computed values of H are much closer to the theoretical ones than the corresponding quantities obtained from the correlation function. Furthermore, part (if not all) of the discrepancy is very likely due to the finiteness of the time sample, not to the Hurst analysis. There is no way that we could possibly use larger r’s in the analysis to check our hypothesis, but the average residence time of the closed state obtained in the numerical simulations for A = 1.3 differs from the theoretical one by some 5%, independently of the number of averages considered. The obvious interpretation is that, due to the finiteness of the sample considered, we are underestimating the tail of the exit time distribution, and this reflects on the numerical value of H. The loss of correlation at larger times should lead to a behavior for R(T) that is more standard when larger T’S are considered, yielding H values that are “closer” to H = 0.5 than the real ones; this is clearly the case. The numerical analysis of the generalized version of the Liebovitch model revealed that the Hurst technique is certainly more efficient than any other possible technique to detect the existence of a power-law behavior of the correlation function C(t) [or a(t)], including the direct evaluation of the correlation function. This is so because the power-law behavior is exhibited in the long-time region where the correlation function is very weak and consequently difficult to distinguish from statistical fluctuations. If the generalized Liebovitch model applied, the Hurst coefficient would provide information on the transition from the short-time region to the long-time region equivalent to those of the square root of the moment (x’(t)). If this held true in general, then it would be more convenient to detect the existence of a power-law behavior by means of the Hurst analysis rather than through the direct evaluation of the second moment, due to the better numerical efficiency. However, we believe that the existence of a power-law behavior

100

ROBERTO

RONCAGLIA

ET AL.

might invalidate the time convolutionless diffusion equation (4.17) and with it the independence of the time scaling of (x”(t)>“” of n. 7.

CONCLUDING

REMARKS

The main results

of this paper

are as follows.

(1) It has been shown that the fractal model proposed by Liebovitch results in an anomalous diffusion if A < 2. (2) It has been proved that the Hurst analysis is a reliable technique to determine “experimentally” whether or not the anomalous diffusion really exists. (3) It has been shown that the power-law decay of the correlation function C(t) might be revealed even in the region A > 2, where standard diffusion is expected to take place in the long-time region. (4) The Liebovitch model has been suitably extended to make it possible for the closed gates to have a finite rather a vanishing statistical weight. There are many problems still open. Note that our numerical analysis is applied to the generalized version of the Liebovitch model rather than to a real ion channel. It is quite encouraging that a recent analysis of real experimental data carried out by Weandle and Liebovitch shows a distinct deviation from standard diffusion with values of H as large as fi = 0.90 [14]. We are convinced that the theoretical results of the present paper can be a useful contribution to the discussion of experimental results of this kind. To take a relevant example, note that the analysis of the generalized version of the Liebovitch model shows that (x2(t)) and R*(t) are virtually identical, thereby suggesting that the generalized Liebovitch model is compatible with the Gaussian diffusion equation (4.17). However, a careful investigation should be carried out on real experimental data to assess whether or not this condition also applies to real systems. We are convinced that the long-time behavior of the correlation function Q,(t) in the long-time regime might accom any the breakdown of the independence of the time scaling of (x”(t))’ P ’ of ~1.This would reinforce the conviction that the dynamics of ion channels has a long memory. Although this paper is merely theoretical, it serves the useful purpose of showing how to conduct a numerical analysis of real experimental data so as to reach the maximum accuracy and to reveal the presence of a power-law behavior in a region that could be mistaken as pertaining to conventional diffusion. We thank the Texas Higher Education Coordinating Board for partial support of this research (Texas Aduanced Research Program, project

FRACTAL

PROPERTIES

101

003594-038) and the European Community for partial financial suppon under contract SCl-CT91-0697 (TSTS). We also warmly thank Dr. Donatella Petracchi for useful discussions about the experimental aspects of the dynamics of ion channels, and Dr. Bruce J. West for introducing us to the subject of fractal dynamics and Hurst analysis. REFERENCES 1 2

3 4 5

6 7 8 9

10 11

12 13 14

M. Barbi and D. Petracchi, Coupled and uncoupled Markov systems: A possible way to distinguish between them, Eur. Biophys. J. 20:345 (1992). D. Colquhoun and A. G. Hawkes, On the stochastic properties of single ion channels, Proc. Roy. Sot. Land. B211:205 (1981). J. Feder, Fracrals, Plenum, New York, 1988. I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, London, 1980. P. Grigolini, The projection approach to the Fokker-Planck equation: applications to the phenomenological stochastic equations with colored noises, in Noise in Nonlinear Dynamical Systems, Vol. 1, F. Moss and P. V. E. McClintock, Eds., Cambridge Univ. Press, Cambridge, 1989, pp. 161-190. S. J. Korn and R. Horn, Statistical discrimination of fractal and Markov models of single channel gating, Biophys. J. 54:871-877 (1988). L. S. Liebovitch, Testing fractal and Markov models of ion channels kinetics, Biophys. J. 55:373-377 (1989). L. S. Liebovitch, Analysis of fractal ion channel gating kinetics: Kinetics rates, energy levels, and activation energies, Math. Biosci. 93:97-115 (1989). L. S. Liebovitch and J. M. Sullivan, Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons, Biophys. J. 52:979 (1987). L. S. Liebovitch, J. Fischbarg, J. P. Koniarek, I. Todorova, and M. Wang, Fractal model of ion-channel kinetics, Biochim. Biophys. Acta 896:173-180 (1987). L. S. Liebovitch, J. Fischbarg, and J. P. Koniarek, Ion channel kinetics: A model based on fractal scaling rather than Markov processes, Math. Biosci. 896:173-180 (19871. B. B. Mandelbrot, The Fractal Geometry of Nature, Freeman, New York, 1983. R. Mannella, B. J. West, and P. Grigolini, A dynamical approach to fractional Brownian motion, Fractals 2:l (1994). S. Weandle and L. S. Liebovitch, Membrane potential fluctuations of human T-lymphocytes have the fractal characteristics of fractional Brownian motion (submitted for publication).