1571
J. Phys. Chem. 1993,97, 1571-1575
Chaos in a Membrane Oscillator Peidong Shen, Jeenok T. Kim, Raima Larter,* and Ken Lipkowitz Department of Chemistry, Indiana University-Purdue University at Indianapolis, 402 N . Blackford Street, Indianapolis, Indiana 46202- 3274 Received:August 3, 1992; In Final Form: October 15, 1992
Aperiodic oscillations in membrane potential in a lipid/alcohol doped Millipore filter are analyzed for the existence of deterministic chaos using standard methods from nonlinear dynamics theory. Previous observations with this system indicated a period-doubling transition into an aperiodic state; both the aperiodic state and the prior periodic states exhibited substantial experimental noise. Here we report calculations of the correlation dimension and Lyapunov exponent for an apparent periodic-chaotic sequence exhibited by this membrane oscillator and find that parameters needed for the implementation of standard algorithms for these two quantities must be chosen carefully for the noisy, unsmoothed data. Nevertheless, reliable values of both the correlation dimension and Lyapunov exponent can be obtained if the dependence on important parameters is treated carefully; the resulting correlation dimension and Lyapunov exponent values are, indeed, indicative of chaos in this system. This work represents the first known example of chaotic behavior in a membrane oscillator.
Introduction
Computer(Mac Ilcx)
In a recent paper' we reported observing oscillations in a membrane constructed by doping a Millipore filter with a 1:l mixture of dioleyl phosphate (DOPH) and oleyl alcohol. The oscillations were observed to go through a sequence of sudden changes, or bifurcations, as the experiment progressed. We concluded that these changes in oscillatory state occurred as a surface emulsion layer became progressively swollen and extended throughout the course of the experiment. The surface layer was determined to be the region of highest electrical resistance and to be critical to the existence of the oscillations. As discussed in more detail in our previous paper, the slow increase in gel thickness can be overcome by making changes in other external parameters over a short time interval during which $he gel thickness is relatively constant. In this paper, we report on our studies of the transitionsbetween oscillatory states in this system. We will be especially concerned with transitions which lead to apparently chaotic behavior. Previous studies2 of DOPH-doped filter systems have been inconclusive about the existence of chaos in these membrane oscillators. With the present results, we can now state definitively that chaos exists in this system. Moreover, the chaotic behavior appears as part of a periodiochaotic sequence, also quite commonly associatedwith chaos. In our previous paper we showed that aperiodic behavior in this system arises by way of a perioddoubling sequence of bifurcations; this prior observation of a well-known route to chaos lends additional credence to our conclusion here.
Experimental Section Membranes were prepared as described in ref 1 and placed into the membrane cell shown schematically in Figure-1. The DOPH used was synthesized by combining oleyl alcohol (Eastman Kodak, purity 60%) with phosphoric acid and refluxing. The solutions bathing the membrane on both sides were 0.1 M KC1. A pressure gradient was imposed by flowing 0.1 M KCl into one sideof the cell; typically, pressures of 30-40 mmHg were imposed acrossthe membrane. A Sage Model 341B syringe pump provided the flow, and the pressure difference was measured with an electronic pressure transducer from Validyne (Northridge,CA). In addition to the imposed pressure gradient, a current of approximately 0.05-2.0 pA (applied with a Keithley current source, Model 224)was also imposed across the membrane. These
-
Recorder
7-1
I Pressure Transducer
I
'
t- I
I
\
Faraday Cage
Figure 1. Apparatus for measuringthe membrane potential. (a) DOPHI oleyl alcohol doped Millipore filter. (b) Ag/AgCl electrodes used to impose the current and measure the resulting potential. (c) Magnetic stirring bars (not used in this experiment). (d) External magnetic stirrer. (e) Solvent outflow.
two external forces triggered the growth of a surface emulsion layer. Once the electrical resistance of the surface layer reached a high enough value, oscillations in the membrane potential could be observed. A Keithley electrometer, Model 617,was used for the latter measurement; data were fed directly to a Mac IIcx computer and analyzed using Labview software (National Instruments).
Theoretical Method Thedata to be analyzed for chaotic behavior consistsof several thousand values of the membrane potential collected over a time interval. In order to analyze this data, we used the method of time delays3 to embed the time series in a two or higher dimensional space. The Grassberger-Procaccia algorithm4 was used to calculate a correlation dimension for each data set. A point Xi in the m-dimensional phase space is given by the vector Xi = (x(i),x(i T ) , ...,x(i ( m- 1)~)) where T is the time lag between consecutive orthogonal axes and x ( i ) is the experimentally
+
+
0022-3654/93/20971571 %04.00/0 0 1993 American Chemical Society
1572 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993
Shen et al.
measured time series data. The correlation function is defined for this embedded data as 4
N-l
N 10
where N is the number of points in the time series, r is a length scale, and 0 is the Heaviside step function. The correlation function C(r) counts the number of pairs of points in the m-dimensional phase space which are less than a distance r apart. For a self-similarattractor C(r)is proportional to rDin a certain range of r values. This power law behavior can be recognized by plotting In C(r)against In (r); a straight line (with slope D)will result if the power law form holds. Alternatively, a plateau in the plot of the slope of In C(r) versus In (r) is also indicative of power law behavior. If the value of the slope at this plateau converges to a constant value with increasing embedding dimension m, the asymptotic value will equal the correlation dimension D,of the attractor. The dynamics of chaotic systems can be further characterized by computing the Lyapunov exponents. An algorithm for calculating the largest Lyapunov exponent A to diagnose chaotic motion for experimental data has been given by Wolf et a1.5This algorithm involves, again, embedding a time series in an m-dimensional space and monitoring thedivergence of trajectories initially close to each other in that phase space. The Lyapunov exponent can then be estimated from
I
0
2
4
2
4
2
4
6
I 0
I
T l m e (second )
where M is the number of replacement steps and d0(rk-,)is the separation distance of two initially close points in the trajectory. After a fixed evolution time At the final distance between the two points evolves to d(rk). Many parameters must be specified to employ the algorithm of Wolf et al., e.g., Ar, the evolution time; the minimum and maximum distance between the two initial points; the allowed orientation error; m, the embedding dimension; and T , the time lag used in the phase space reconstruction. We have carried out an analysis of the dependence of the Lyapunov exponent on these parameters and have found that the evolution time At is especially important in diagnosing chaotic behavior in this membrane oscillator and, possibly, in other physical systems. This point will be further clarified in the discussion.
Results and Discussion Approximately 30 min after the pressure gradient is applied to the membrane, a thick, opaque DOPH/oleyi alcohol/water (D/O/W) emulsion appears on the low-pressure surface of the membrane. This gel-like layer grows as the experiment proceeds; the growth of this surface layer triggers a dramatic rise in the electrical resistance of the membrane. When the electrical resistance reaches several MQ, imposition of a small current (0.052.0 PA) results in oscillations in the membrane potential. The growth of the gel also causes a buildup in the pressure gradient across the membrane; oscillations are seen when this pressure gradient is approximately 30-40 mmHg. The formation of the high electrical resistance gel-like layer is the key event which triggers oscillations in the membrane potential. The specificoscillatory state observed depends on the current (I),the pressuregradient (AP), and theextent of expansion of the gel layer. As the experiment proceeds, the gel layer continues to expand and bifurcations in oscillatory state are observed. As described in ref 1, it was found that if the current I is held constant but the gel layer continued to expand, the oscillatory state would undergo transitions from simple periodic behavior to complex periodic to an irregular, possibly chaotic,
Figure 2. Time series of membrane potential measurements at I = 0.5 p A and AP = 3 0 4 0 mmHg. The membrane was made by using a 1:l (by volume) mixture of DOPH and oleyl alcohol. The amount of D/O adsorbed to the membrane was typically 6-9 mg/cm2. The solutions bathing themembraneonbothsides were0.1 M KCI. Different oscillatory states a, b, c, and d were recorded as the gel expanded.
state. Figure 2 shows the oscillatory time series for a typical experiment in which I = 0.5 MA. Here the thickness of the gel layer is the bifurcation parameter which causes changes in the oscillation pattern from the periodic state a, to a second state b, then to the aperiodic state c, and returning to an apparently periodic state d. Sudden changes in oscillatory state without deliberate alteration of an externally controllable parameter have also been observed in the oscillatoryelectrodissolutionof copper.6 These sudden changes were attributed to the growth of an oxide film on the copper surface; the gel layer here is analogous to the oxide film, although the specific chemical mechanism is clearly different in the two cases. The data shown in Figure 2 have not been smoothed or averaged in any way and, hence, contain a fair amount of noise typical of these membrane oscillators. The membrane potential readings were saved directly to computer disk at intervalsof 0.004s. Figure 3 shows the corresponding two-dimensional (V(t),V(r + T ) ) phase portraits constructed by the method of time delays with T = 0.008 s. Although the data are quite noisy, it is possible to see what appears to be a period-doubling bifurcation when going from 3a to 3b. The second oscillatory state is so noisy, however, that it is impossible to determine from the phase portrait alone whether it is periodic or not. The data in 3c show what appears to be a fully-developed strange attractor. Beyond this possible chaotic state, a more ordered, possibly periodic, state appeared as shown in 3d. Becauseof the high amount of noisein thesedata it is imperative to carry out correlation dimension and Lyapunov exponent calculations in order to ascertain with certainty whether thestate 3c is truly chaotic. Figure 4 shows the results from the correlation dimension calculation as a plot of In C(r)versus In (r) for the four states in Figure 2. The embedding dimension was varied over the
The Journal of Physical Chemistry, Vol. 97, No. 8, 1993 1513
Chaos in a Membrane Oscillator
-: z
2
4
I
0
- 5
- 3
- 1
- 5
- 3
- 1
10
c
a
I
+
a
I
10
12
10 -I
I
I"
I
1
.a
I
I
-
0 - 4
- 3
- 2
- 1
Figure 5. The corresponding slope of In C(r) versus In (r) in Figure 4. In the middle range of r in each plot, a 'plateau" can be seen (indicated by arrow). A blowup of this plateau region is shown in Figure 6 .
..-...-._..
P)
a
I 0.6 -3.8
-3.3
-3.0
6
7
9 '0
Figure 6. The slope of In C(r)versus In (r) for the blowup of the plateau regions in Figure 5 . The slope of all four states is seen to converge with increasing embedding dimension. Figure 4. Correlation dimension calculations from eq 1. Plots of In C(r) versus In (r) for the four states a, b, c, and d. The number of data points used in these calculations was N = 3000. The numbers labeling the curves in the plot indicate the embedding dimension, chosen to cover the range 2-1 I .
range m = 2-1 1, but it is not possible to determine from these plots whether the correlation dimension has converged. Plots of the slope of the graphs in Figure 4 are shown in Figure 5. These results show a great deal of scatter at both low values of r and high values of r; however, in the middle range of r values, a 'plateau" can be seen (indicated by the arrow). A blowup of this plateau region is shown in Figure 6. In all four cases, the value of the slope in this midregion is seen to converge with increasing embedding dimension. Thegeneral dependenceof the correlation
dimension on these parameters is similar to that reported by Harding and Ross7 for a different system. Figure 7 shows the convergence of the correlation dimension to the asymptotic value for the four sequential states. In each case, increasing the embedding dimension m beyond the last value shown did not lead to any further changes in the asymptote. The asymptotic value is very close to 1 for the initial state, identified tentatively as periodic (a). It increases to a value between 2 and 3 for the intervening aperiodic states (b and c) and then returns back to a value near 1 for the next periodic state (d). The relatively low value of the correlation dimension for the middle two states is indicative of chaos just past the transition. This periodicchaotic sequence is typical of the results seen in other nonlinear oscillators.
Shen et al.
1574 The Journal of Physical Chemistry, Vol. 97, No. 8, 1993 1.2
2.3
a)
0.5 0.4 0.6
1.o
L 5
3
2.0 2
4
6
m
0.3
-
OS
-
8 1 0 1 2
m
6.P
Figure 9. Lyapunov exponents calculated from eq 2 versus evolution time At for the Henon map of Figure 8. For low evolution times At (Le., the first 2 points), A = 0.60 agrees very well with the accepted value8of
a n 1.8
0.603.
-
-
0.8
1.4 2
4
6
8
10
12
0.4
I 2
6
4
8
10
12
m m
Figure 7. The values of slope in Figure 6 versus embedding dimension. The convergenceof the correlation dimension to the asymptotic value can be seen for the four sequentialstates. Theconvergedcorrelation dimension D, is 1.2 in a, 2.3 in b, 2.6 in c, and 1.5 in d.
-- .
0.04
0.00
1-
X
0-
-1
-
L
-2 0
0.08
0.12
Evolution t h e (second)
50
100
T h e ( second )
Figure 8. The chaotic time series for the Henon map with a = 1.4, and 6 = 0.3. The time interval (between points n and n 1) was taken as 1 s here. The Henon map is a pair of nonlinear equations: X,+l = Y, 1 - aXn2 and Y,+, = 6Xn. Once initial values X, and Yo are given, a 'time series," Le., values of X,,and Y,,,can be obtained by iterating the map.
+
+
The calculations of the Lyapunov exponent X were found to be relatively independent of most parameters used in the calculations but to depend heavily on the evolution time At. It was difficult to find a definite range of At values in which the calculated value of X was constant. To test this result further, a simple but well-understood model, the Henon map, was chosen to study the dependence of the Lyapunov exponent calculation on the choice of parameters. Because of the discrete nature of our experimental data a system of difference equations (such as the Henon map) provides the best means of comparison. The chaotic time series for the Henon map is shown in Figure 8. The accepted value of the Lyapunov exponent* is X = 0.603 (where Xis defined in terms of d(t) = do2*'). The units of time, i.e., the value off between the nth point and the ( n 1)th point, were taken as 1 s. The calculated value of X for Figure 8 is plotted as a function of the parameter At, the evolution time, in Figure 9. At small evolution times, the calculated X agrees with the accepted values of h very well. The Lyapunov exponent was also calculated for our time series a, b, c, and d of Figure 2. The resulting values of A as a function of evolution time are shown in Figure 10. Reliable values of X
+
Figure 10. Lyapunov exponents calculated from eq 2 versus evolution time At for the four experimental time series. Reliable values of A were obtained at low evolution times (first 2 or 3 data points). The states a and d have A values which are negative and close to zero, which indicates periodicstates. Thestates bandc have positivevalueof A,whichindicates chaotic states. In addition, A increases from b to c as the system becomes more highly chaotic.
TABLE I: Lynpuaov Exponents Calculated from the Wolf Algorithm for the DOPH Membranes time series
A
time series
A
a b
-0.02 0.08
C
0.25 -0.007
d
are obtained at low evolution time (which corresponds to the first 2-3 data points). This is consistent with the dependence on evolution time seen for the Henon map. The states a and d have A values which are negative and close to zero, which indicates they are periodic states. The states b and c have positive values of X which indicates chaotic states. In addition, the computed value of X increases from b to c as the system becomes more highly chaotic. The time interval between data points in the experimental data was 0.004 s. Renormalizing this interval so it is the same as the Henon map (i.e., 1 s instead of 0.004 s) we find that the absolute values of the computed Lyapunov exponents become quite reasonable for this low-dimensional chaos. The values of the renormalized Lyapunov exponents at small evolution times are given in Table I. The values of and h, are typical for a chaotic system just past the transition to chaos.
Conclusions Previous studies of similar doped-membrane systems have led to inconclusive evidence that the aperiodic behavior observed in these systems is an example of chaos. Hayashi, et al.2 showed several examplesof aperiodicoscillations observed with a Millipore filter doped only with DOPH and an unknown amount of oleyl alcohol. A calculation of the correlation dimension revealed values
Chaos in a Membrane Oscillator of 5 or more for these latter data. Correlation dimensions as large as these may be an indication that the aperiodic states were actually just random experimental noise; chaos cannot be definitively ruled out for their results, however. We calculated the correlation dimensions and Lyapunov exponents for our DOPH membrane system via the algorithms of Grassberger-Procaccia and Wolf et al. without subjecting the experimentaldata to presmoothing or averaging. We found that the Lyapunov exponent depended heavily on the evolution time, a parameter needed for implementation of the Wolf algorithm. The Henon map was used to determine the best range of values for this parameter and others. Reliable Lyapunov exponents could be obtained for our data using low values of the evolution time. The resulting values of 0.08 and 0.25 for the aperiodic states are typical for a chaotic system. The values of the correlation dimension are somewhat higher than would be expected; for example, for the state in Figure 2a a correlation dimension slightly greater than 1 was found. A truly periodic state would be expected to have a correlation dimensionof exactly unity. It has been showng that noise tends to impart a slight increase to the computed values of the correlation dimension; this is a possible explanation for our results. A pssible source of noise in this system is spatial nonuniformitiesor slight variations in thickness of the gel, although these were not detectable in this experiment. It would not be surprising if spatial patterns were to exist in this type of system.
The Journal of Physical Chemistry, Vol. 97, No. 8, 1993 1575 The values of 2-3 found here for the aperiodic oscillations are most likely to be observed just past the transition to chaos. We previously observed' what appeared to be a period-doubling sequence of states leading up to an apparently chaotic state. This work provides definitive evidence that this latter state actually is deterministic chaos, and not just an example of experimental noise in the data, and furthermore, constitutes the first proof of chaotic behavior in a membrane system.
AcLnawledgwnt. Theauthors gratefullyacknowledge support of this work by the donors of the Petroleum Research Fund, administered by the American Chemical Society,and the National Science Foundation under grant CHE-8913895. We also thank Professor K. Toko of Kyushu University for helpful discussions.
References rad Notes T.;Larter, R. J . f h y s . Chem. 1991. 95, 7948. (2) Hayashi, K.; Toko,K.; Yamafuji, K. Jpn. J . Appl. f h y s . 1989,28, 1507. (3) Packard.N.;Crutchfield, J.P.;Farmer, J. D.;Shaw,R.S.fhys.Reu. Left. 1980, 45, 712. (4) Grassberger, P.; Procaccia, I. fhysicu 1983, 9 0 , 189. (5) Wolf, A.; Swift, J. B.; Swinney, H. L.; Vastano, J. A. fhysicu 1985, lbD, 285. (6) Bassett, M. R.; Hudson, J. L. J . f h y s . Chem. 1988, 92,6963. (7) Herding, R. H.; Ross, J. J . Chem. f h y s . 1988, 89, 4743. (8) Simo, C. J . Sfuf. f h y s . 1979, 21, 465. (9) Geest, T.;Olsen, L. F.; Schaffer, W. M. In preparation for Biochim. Biophys. Acru, 1992. (1) Kim, J.
