J. Phys. Chem. 1992, 96, 4931-4934
4931
Propagation Failure in Arrays of Coupled Bistable Chemical Reactors Jean-Pierre Laplante Department of Chemistry and Chemical Engineering, Royal Military College of Canada, Kingston, Ontario, Canada K7K 5 u )
and Thomas Emeux* Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, Illinois 60208 (Received: January 27, 1992)
We consider a onedimensional array of coupled stirred tank reactors. Each reactor operates with the bistable chlorite-iodide reaction. The reaction admits two stable steady states (A and B). At time zero, all reactors are at steady state A except the first reactor which is at steady state B. Provided the coupling (exchange rate) between each reactor is sufficiently large, a propagating wave joining A to B is initiated. We determine the speed of the front and show that it fails to propagate if the exchange rate is below a nonzero value. In addition, we investigate the time history of the front and identify three distinct stages. First, we note a small time regime dominated by a diffusion mechanism; second, we note a front propagating with almost constant speed; third, we note a boundary layer regime corresponding to the end of the front and the approach to the final uniform steady state. We explain all three stages by analyzing a simple model.
1. Introduction Evidence of propagation failure as well as other conduction irregularities (delay, partial reflection) has been found experimentally with cardiac tissues and other excitable membrane systems. Propagation with a nonconstant velocity may result either from a weak coupling between the cells or a nonuniform distribution of the parameters in the population. Conduction anomalies may lead to fatal diseases, and computer simulations of nerve conduction have modeled various mechanisms leading to propagation Detailed experimental observations of this phenomenon are technically difficult with biological tissues. In this paper, we investigate a system of coupled bistable chemical reactors as a simple model of a chain of cells. We propose to determine how a front joining the two stable steady states may fail to propagate. Coupled stirred flow reactors (CSTRs) have previously been investigated for their instabilities. The experimental system usually corresponds to assemblies of two4+ to seven r e a c t ~ r s . ~In particular Marek and c o - w ~ r k e r sand ~ ~ more ~ recently Marmillot et al.7dhave found a large variety of symmetrical and asymmetrical stationary structures suggesting the possibility of cascading bifurcations. In addition, bifurcation of periodic solutions might o c c ~ r . In ~ this ~ , paper, ~ ~ we concentrate on the case where the reaction taking place in the reactors is bistable. We use the chlorite-iodide reaction which has been extensively studied by Epstein and co-workers (for a review, see ref 8). In a CSTR, this reaction is known to exhibit multiple steady states over a wide range of parameters. By contrast to previous studies on instabilities appearing in coupled C S T R S , ~our ~ main interest is to determine the conditions for a propagating wave. In a continuous bistable reaction-diffusion system, a wave propagating well away from the boundaries can be characterized by its constant speed c which depends on the kinetic parameters and the diffusion coefficients. In a system of coupled bistable reactors, coupling only appears at the boundaries of each reactor and a wave may propagate with a relatively small velocity. Our experiments clearly indicate that the speed of the front is continuously changing with time. We have found three different stages corresponding to the small time, intermediate, and final regimes of the propagating front. During the intermediate regime, the speed is almost constant. As the exchange rate decreases and approaches a nonzero value, we have found that its value approaches zero. No traveling waves are observed below this critical value. The plan of the paper is as follows. Section 2 describes the experimental results. Section 3 is devoted to the study of a simple model. We show that the computer simulation of the propagating wave reproduces the main experimental observations. We discuss our results in section 4.
2. Experiments A. Reactor Assembly. As illustrated in Figure 1, we consider 16 reactors connected in a linear geometry. The coupling between each reactor is by mass exchange. The reactors were machined out of acrylic and within our experimental accuracy considered as identical. Their internal volume was measured to be 5.80 f 0.05 mL. In our experiments, each of the reactors was fed at a constant flow rate with the chlorite and iodide reagents, through two 16channel peristaltic pumps (Ismatec). To avoid any artifact due to mixing, the reactants were premixed in a small volume T-union just before entering the reactor^.^^^ The standard deviation of the flow rate, from tube to tube, was typically found to be of the order of 2-3%. A slow drift in the average flow rate was also observed as the tubes were aging. The tubes were therefore periodically recalibrated or replaced to improve reproducibility. The stirring was magnetic and, through a chain and sprocket system, identical in all reactors. The experiments reported in this paper were all carried out at 850 rpm and at room temperature. In order to implement mass exchange between adjacent reactors, two more 16-channel peristaltic pumps were used. As seen in Figure 2, the reactors were split in two groups, the odd and the even reactors. These two groups were gathered as close as passible, on opposite ends of the peristaltic pumps. In order to minimize the length of the exchange path, the exchange pumps were positioned above and below the reactors. One pump provides the mass transfer in one direction whereas the other takes care of the transfer in the opposite direction. From connection of the in and out of each of the flow tubes to the appropriate reactors, the linear geometry shown in Figure 1 was realized. Again here, special attention was paid to the calibration of the exchange tubes to ensure that the flow rate was identical in both directions. Note that this method had the disadvantage of introducing an inherent delay in the transfer of material between adjacent reactors. This delay is equal to the time required for the solution to transit from a given reactor to its neighbors and is inversely proportional to the exchange rate. In the experiments described herein, this delay was kept below 5 s by using tubes of various inside diameters. B. R e d & In order to study a propagating wave in this system of coupled reactors, we have used the following procedure. With the exchange reactors between adjacent reactors turned off, reactors 2-16 were identically prepared in a steady state on the lower branch of the chlorite-iodide butability curve (see Figure 1). The essential characteristics of the wave propagation were found to be relatively insensitive to the exact location of this steady state, provided it remains located within the hysteresis limits. On the other hand, reactor 1 was used as a constant concentration boundary condition and continuously fed with a concentrated
0022-3654/92/2096-493 1%03.00/0 0 1992 American Chemical Society
4932 The Journal of Physical Chemistry, Vol. 96, No. 12, 1992
-
16
COUPLED REACTORS
14
More "'Reactors "'
highI [chlorite] (constant concentration)
Laplante and Erneux
12
I no flux
10 8
6 4
0
500
[170 Figure 1. Schematic representation of the experiment. Reactors 2-1 6 are prepared in a steady state on the lower branch of the chlorite-iodide bistability curve whereas reactor 1 is fed with a chlorite solution only. As the exchange between reactors is turned on, a wave front corresponding to the transition of reactors 2-16 from the lower to the upper branch of steady states is initiated. On the bistability diagram [I-] refers to the steady-state concentration of iodide in the reaction, whereas [I-],, refers to the concentration of iodide resulting from the mixed reactant feed.
Top View
2000
(9)
Figure 3. Wave propagation through a linear array of bistable chlorite-iodide media. The propagation of the wave front is monitored through the color change (blue to colorless) that accompanies the transition from the lower to the upper branch of steady states of the chlorite-iodide reaction (see Figure 1). The time of change refers to the time at which the color change takes place in a given reactor, time 0 being defined as the time at which the exchange was turned on. For the experiments, reactors 2-15 were all fed at identical flow rates, with the M in pH 3.5 buffer following solutions: (1) NAI, 1.666 X (Na2S04/H2S04);(2) NaC102, 4.500 X 10" M in 0.001 M NaOH. A 0.25 g/L amount of soluble starch was added to the NaI feed. The residence time (in absence of mass exchange) was kept constant at 170 s. Reactor 1 was fed with a 0.02 M solution of NaC102 in 0.001 M NaOH. Note that in order to prevent contamination of the chlorite solution from the exchange coming from reactor 2, the exchange flow from reactor 2 was sent to waste. I
Even Reactors
Odd Reactors
1500
1000
Time of Change
0.012
-
0.002
-
I
I
I
I
I
I
I
6
7
8
9
1
\
Odd Reactors
16 Channel
Peristaltic Pumps
Even Reactors
Side View
Figure 2. Physical configuration of the CSTR assembly. Reactors are split in two groups. Sixteen-channel peristaltic pumps are used to implement mass transfer according to the linear geometry shown in Figure 1.
solution of sodium chlorite reagent. In order to prevent any contamination of this boundary solution from the exchange flow of reactor 2, the exchange flow from 2 to 1 was sent to waste. Note that this choice of boundary condition corresponds to maintaining reactor 1 at a fmed steady state (Le., [I-I0 = 0) on the upper branch of the bistability curve (see Figure 1). The steady states corresponding to the upper and lower branches of the bistability curve can easily be distinguished by adding soluble starch to the reactants in flow. On the lower branch, the solution takes a deep blue color characteristic of the presence of a high concentration of the triiodide ion. On the upper branch, the solution is almost colorless.* A propagating wave was initiated by turning on the exchange flow between the reactors. Even though each of the reactors was fitted with a redox electrode, it was found more convenient to monitor the propagation of the wave by using the color change that accompanies the transition from the lower to the upper branches of steady states. Figure 3 illustrates the results of a group of four experiments carried out at different exchange rates. These results clearly show that the speed of the wave front is continuously changing with time. Three different stages can, however, be identified: a small time (diffusion-like) regime, an intermediate phase, and a final regime during which the front apparently accelerates. We investigate these three regimes in the next section by using a simple model. We note that the intermediate regime can be characterized by the fact that the wave-front velocity is
d
\ I
0.000
4
5
10
Exchange r a t e (mL/min) Figure 4. Propagation failure in a linear array of coupled bistable reactors. The stationary wave-front velocity is clearly seen to approach zero as the exchange rate approaches the critical value d*. These velocities were obtained from the slope of the curves shown in Figure 3 at the inflection point. Experimental conditions are the same as in Figure 3. The solid line corresponds to expression 12 where d* and p have been obtained so that (12) is the best bit of the experimental data. We obtain d* = 6.180 mL/min and p = 0.01498 reactor/[s (mL/mir~)'/~]. The dashed line corresponds to the parabolic law (9) with d* = 6.259 mL/min and y = 0.005 871 reactor/[s (mL/min)1/2]. These values provide the best fit. Note that the range of exchange rates used in the results shown above corresponds to exchange times of the order of 45-60 s, which is about 3 times smaller than the residence time in the absence of exchange. Propagation failure was also observed in other series of experiments where the typical exchange time was significantlysmaller than that given above.
almost constant. The propagation in this regime is then similar to what would be observed in an infinite array of discrete cells. Of particular physical interest is how this velocity depends on the exchange rate (see Figure 4). We have determined the velocity from the slope of the curve shown in Figure 3 at the infection point. From Figure 4, we note that the velocity decreases with the exchange rate and becomes zero at a critical (nonzero) value. Below this critical value, the front is stationary and separates the
The Journal of Physical Chemistry, Vol. 96, No. 12, 1992 4933
Arrays of Coupled Bistable Chemical Reactors reactors that have made the transition to the upper branch from those which, if unperturbed, will remain on the lower branch of the steady states. This failure of the wave front to propagate through the entire array of reactors is called propagation failure. As shown mathematically by Keener,lo propagation failure is a property of the discrete system and cannot be observed in a continuous, one variable, reactiondiffusion system if the parameters are uniformly distributed. If a parameter is not constant throughout the medium, propagation failure is possible in a continuous reactiondiffusion system and plays an important role in population genetics or neurobiology.'JL 3. Theory In this section, we propose a simple model for coupled bistable reactors and investigate the three stages of the moving front observed experimentally. These three stages are (i) a small time regime dominated by a diffusion-like mechanism, (ii) an intermediate regime characterized by a almost constant speed of propagation, and (iii) a final approach to the uniform steady state where the front accelerates. Our model corresponds to coupled Nagumo equations and is the simplest model for a system of coupled bistable CSTR. Because of its simplicity, it has been analyzed in detail both numerically and analytically.L0J2 In this section, we show that the main features of the front evolution can be explained by this model. The model consists of the following equations: du,/dt = d(u,, - 2u, + u,+~) +flu,), 1 I n I N (1)
duN/dt d(UN-1 - uN) + f l ~ N ) where the nonlinear function fiu) is defined by flu,a) = u(u - l)(a - u) 0I a I 1
(2)
(3) The model depends on two parameters (a and d). d is the control parameter and corresponds to the exchange rate. Equations 1-3 admit three uniform steady states given by u, = 0, u, = a, and u, = 1. We now propose to determine a traveling wave solution of eq 1-3. Note that f'(0) < 0 and f'(1) < 0. Thus, the uniform steady-state solutions u, = 0 and u, = 1 (1 I n I N) are stable solutions of eqs 1-3. We are interested in finding a traveling wave solution which is initiated by the following boundary and initial conditions: uo=1
t 1 0
(4)
l l n l N and t = O (5) Provided that d is sufficiently large, a numerical study of eq 1-5 shows that u,(t) is changing from u = 0 to u = 1 at a fixed time t = t,. We define t , as the time corresponding to u, = u* where u* is an arbitrary reference value (0 < u* < 1). In Figure 5, we represent t , as a function of n (solid line). The values of the parameters are a = 0.2, d = 1, and u* = 0.1. The curve summarizes the time history of the propagating wave. As for the experimental curves shown in Figure 3, we observe three different stages: (i) the small time solution, (ii) the almost constant speed solution, and (iii) the boundary layer solution as the wave approaches a uniform steady state u, = 1. A. Small Time Solution. The small time solution is dominated by a diffusion mechanism. In the Appendix, we construct the small time solution and show that the leading approximation for t , is given by
l57
i
8
,, , ,, /
"'
2.32 (d tn)
: L
0
2
-1
0 0
10
20
30
40
time of change t, Figure 5. Time history of the moving front. The curve has been obtained numerically from eqs 1-5 and represents the reactor number n as a function of the time of change f , (0 9 n IN). This time is defined by the condition I&,,) = u s . The values of the parameters are a = 0.2, d = 1, u* = 0.1, and N = 15. All quantities are dimensionlessquantities. The dashed line corresponds to the small time approximation given by (7).
B. Constant Speed Solution. Away from the boundaries, the propagating wave behaves essentially as a function of a single variable, the traveling coordinate, z = t - n/c, where the speed c is a constant determined by the parameters d and a. KeenerLo analyzed this problem and showed that
-
c
= ( K - a)(24'/2
(8) as d -. In addition, KeenerLodemonstrated that the speed c approaches zero as d approaches a constant d+. Using bifurcation arguments which are based on the existence of a limit point of a nonuniform steady state, it can be shown that
-
c = y(d
-
- d*)1/2
(9)
as Id - d*l 0 (y is a fixed number). In ref 12, Erneux and Nicolis analyze the limit a 0, d = O(a2) and N fixed. They obtain
u,=O
tn = n2/(4C2d)
(6)
C represents a constant which is the root of the equation u* = erfc (C) where erfc ( x ) is referred to as the error-function com-
~ 1 e m e n t . I ~If u* = 0.1, (6) becomes f ,= 0.186(n2d1) (7) This function is represented in Figure 5 by a dashed line. For progressively higher values of d, the agreement between the solution (7) and the exact numerical solution is over a longer period of time.
-
where d* = 1/4a2.This function admits asymptotic limits which are in qualitative agreement with the limits (8) and (9) (c 2 ~ - ' d ' /as ~d and c 7r-l(d- d*)I/*as d d*, respectively). Expressions 9 and 10 suggest two approximative curves for the experimental data. The dashed line in Figure 4 corresponds to the parabolic law (9) with parameters d* and y. The solid line corresponds to the expression
-. -
-
with a = 2(d*)1/2and parameters p and d*. Expression 11 is suggested by (10). For both curves, the two parameters are determined numerically to give the best fit with the experimental data. C. Final Approach to the Steady State = 1. The final stage of the front is the approach to the uniform steady state. It corresponds to a complex boundary layer solution where the amplitude of the front progressively decreases (see Figure 6). It is the decrease of the amplitude of the front which explains its apparent acceleration. 4. Discussion We have shown that the three stages of the evolution of the front which have been observed experimentally can be described
4934
The Journal of Physical Chemistry, Vol. 96, No. 12, 1992
Laplante and Erneux U,+I(f)
u,
"2
= u(i$+7t-I/2,t)
U,l(t) = We now investigate the limit T (A3) and (A4), we obtain
1 .o
-
(A3)
(A41 0, with t small but fmed. From
U(~Tt-'/2,t)
0.8
u,+I(t) = u ( [ , f )
1 + T f - l / ' t + ( [ , f ) + -T2t-'UIt(f,t) + ... 2
0.6
u,l(t) = u ( [ , t )
- Tt-I/*Ut(t,t)
0.4
1 + -T*t-'u,e([,t) + ... 2
Then, with (A2), (AS), and (A6) and with
T
= &l/2 eqs 1-5 lead to the following problem as T 1 1 UI - t t u t =flu) + ,Ut{
-
7
0.2
0.0
u(0,t) = 1, ut(m,t) = 0, and
n Figure 6. Moving front. The solution u,(t) of e q s 1-5 is represented as a function of the reactor number n for different times. Note that the amplitude of the front decreases near the second boundary of the system ( t = 25,30, and 35). The values of the parameters are the same as those in Figure 5.
by a simple model of coupled bistable reactors. Thus,the evolution of the wave does not depend on the particular details of the chlorite-iodide reaction but rather the bistability property of the reaction and the diffusion-like coupling mechanism. Approximate solutions were obtained for the small time solution and the constant speed solution. The small time solution shows that the behavior of the front does not depend on the reaction in first approximation and is useful if one is only interested in the effect of the coupling mechanism. The constant speed solution is also useful because it does not depend on the boundary conditions in first approximation and simple approximations for the speed can be obtained either analytically or numerically. The dispersion of the experimental results for the speed as a function of the exchange rate (Figure 4) is unfortunately too large to assess which of the expressions (9) or (12) provides the best description in this case. Note that this dispersion depends not only on the accuracy of all the parameters (such as the pump rates of each reactor) but as well on the pronounced stochasticity of the chlorite-iodide r e a ~ t i 0 n . l ~ Acknowledgment. We thank G. Dewel and A. Arneodo for stimulating discussions. J.-P.L. thanks the Ministry of National Defence of Canada for supporting this project (Grant. No. FUHDS). The work of T.E. was supported by the US.Air Force Office of Scientific Research under Grant No. AFOSR-90-0139 and the National Science Foundation under Grant No. DMS900 1402.
Appendix. Small Time Solution The basic idea for the small time solution is that pure diffusion is the main mechanism. This suggests that u,(t) = u(t,t) is a function o f t and t as f 0 where [ is the diffusion similarity variable defined a s
-
[=
nTt-l/2
(AI)
In (AI), T is a constant which will be specified later. Definition A1 implies the chain rule dun -= dt
1
ztue
as well as the following expressions for u,+l(t) and u,,(t)
('42)
(A5) (A6)
defined as (A7)
0:
u(m,O) = 0
(A81 (A9)
We solve eqs A8 and A9 by seeking a small time solution of the form u = u(t,t) = u&) + tu1(() + ... (A101 Substituting (A10) into (A8) and (A9), we obtain in first approximation the following equation for uo(t): -I/zEuot
uo(0) = 1, u,(..)
= Uoet = uot(m) = 0
(All) ('412)
This equation is an ordinary differential equation, and two successive integrations lead to the solution uo = erfc (5/2) (A13) where erfc (x) is referred to as the error-function ~mmp1ement.l~ In terms of the original variable, we have found that
--
u,
as t
= erfc (n/(2d1/211/2))
-
(A141
0 and for d sufficiently large (since from (A7) the limit T 0 implies the limit d a). The critical time t = t , corresponding to u = u* satisfies the relation u* = erfc (n/(2dlk,1/2)) (A1 5 )
Equivalently, we may write t , as a function of n t , = n2/(4Qd)
('416)
where C is a constant satisfying u* = erfc (C). Registry No. CIOz-, 14998-27-7; I-, 20461-54-5.
References and Notes (1) Keener, J. P. J . Theor. Biol. 1991, 148,4942. (2) Bell, J. In Reaction-Diffusion Equations; Brown, K. J., Lacey, A. A., Eds.; Claredon Press: Oxford, U.K., 1990; p 95. (3) Rinzel, J. Ann. N.Y. Acad. Sci. 1990, 591, 51. (4) Boukalouch, M.; Elezgary, J.; Arneodo, A.; Boisonade, J.; De Kepper, P. J. Phys. Chem. 1987, 91, 5843. ( 5 ) Bar-Eli, K.; Reuveni, S . J . Phys. Chem. 1985,89, 1329. (6) (a) Crowley, M. F.; Field, R. J. J . Phys. Chem. 1986, 90, 1907. (b) Crowley, M. F.; Epstein, I. R. J . Phys. Chem. 1989, 93, 2496. (7) (a) Stuchi, I.; Marek, M. J . Chem. Phys. 1982,77,2956. (b) Dolnik, M.; Marek, M. J . Phys. Chem. 1988, 92, 2452. (c) Dolnik, M.; Padusaka, E.; Marek, M. J. Phys. Chem. 1987, 91,4407. (d) Marmillot, P.;Kaufman, M.; Hervagault, J. F. J . Chem. Phys. 1991, 95, 1206. (8) De Kepper, P.; Boissonade, J.; Epstein, I. R. J. Phys. Chem. 1990, 94, 6525.
(9) Menzinger, M.; Giraudi, A. J . Phys. Chem. 1987, 91, 4391. (10) Keener, J. P. SIAM J. Appl. Math. 1987, 47, 556-572. ( 1 1) (a) Rinzel, J. Lecr. Appl. Math. ( A M S ) 1981, 19, 281. (b) Pauwelussen, J. J . Math. Biol. 1982, 15, 151. (12) Erneux, T.; Nicolis, G. Submitted for publication. (13) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions; Dover: New York, 1972. (14) Epstein; Nagypal J . Chem. Phys. 1988, 89, 6925.
