Realistic Parameters for Simple Models of the Belousov–Zhabotinsky

ARTICLE pubs.acs.org/JPCA

Realistic Parameters for Simple Models of the Belousov
Zhabotinsky Reaction Jerzy Gorecki,*,†,‡ Jan Szymanski,† and Joanna Natalia Gorecka§ †

Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland Faculty of Mathematics and Natural Sciences, Cardinal Stefan Wyszynski University, Dewajtis 5, Warsaw, Poland § Institute of Physics, Polish Academy of Sciences, Al. Lotnikow 36/42, 02-668 Warsaw, Poland ‡

ABSTRACT: We use experimental results to estimate the values of parameters of simple models describing the time evolution of the Belousov
Zhabotinsky reaction proceeding in droplets surrounded by hydrocarbons. The equations with fitted parameters correctly describe the period of oscillations for a large class of experimental conditions at which the reaction is performed.

1. INTRODUCTION An increasing interest in applications of chemistry to information processing has been noticed in recent years.1
4 This trend can be explained with quite obvious observation that chemical reactions are responsible for efficient information processing tasks performed by living organisms. Investigation of chemical information processing can help in better understanding such phenomena. It is also expected that studies on biological and chemical computing can have an inspiring role and lead to new algorithms for complex tasks like optimization5 or image recognition.6 A typical nerve system is built of nonlinear elements (neurons) that are linked together and form a network. An excitation of one neuron can influence all of those with which it is linked. Experimental studies on real neurons are difficult and costly.7 Therefore, it is not surprising that chemists have recently focused attention on their simpler substitute, droplets containing reagents of a reaction that exhibits complex, nonlinear dynamics separated by surrounding, immiscible liquid. Droplets with a water solution of reagents of the Belousov
Zhabotinsky (BZ) reaction surrounded by hydrocarbons are the most common version of such a system.8 However, droplets in hydrocarbons are quite unstable and merge easily. In order to increase droplet stability, the hydrocarbon solution of surfactants8,9 or lipids10 is used. In such a case, lipids form a monolayer at the boundary between water and hydrocarbons. The presence of a monolayer stabilizes the droplets mechanically, and their structures remain unchanged for hours. For information processing applications, it is important that droplets can interact via exchange of reagents through the monolayer and the surrounding hydrocarbons. These features combined with a long time within which the medium inside of droplets remains active make a BZ droplet system an interesting candidate to study information processing phenomena. It has been shown4 that in information processing applications, the structure of the medium defined as the geometrical r 2011 American Chemical Society

distribution of regions characterized by different chemical dynamics plays an equally important role as the properties of the chemical reactions involved. In the system of BZ droplets, information can be coded in the presence of excitations at selected droplets or in the spatiotemporal correlations between excitations. For example, by arranging droplets in a specific manner, one can obtain structures operating as the basic logic gates.11,12 Theoretical studies on structures interesting for information processing operations are supported by numerical simulations based on models of chemical dynamics. Simple reaction models like Oregonator13 are widely used. However, usually, the parameters of such models are selected quite arbitrarly. The required character of time evolution predicted by the model is used as the main argument to justify the particular choice of its parameters. Therefore, it is quite hard to translate the parameters of such models into real experimental conditions. In this paper, we describe two still simple but yet realistic models that can be used to simulate the time evolution of chemical processes occurring in droplets with the ferroin catalyzed BZ reaction. We fix their parameters to obtain the best agreement with a large class of experimental results. The proposed models take the initial concentrations of reagents as input parameters; therefore, the results of simulations can be directly used to predict experimental conditions in which the calculated information processing operations should appear.

2. SIMPLE MODELS FOR OSCILLATORY EVOLUTION Our experiments have shown that droplets containing a solution of reagents of the Belousov
Zhabotinsky reaction surrounded by a solution of phospholipids in hydrocarbons are stable considering both chemical and mechanical criteria.10 The Received: April 7, 2011 Revised: July 1, 2011 Published: July 12, 2011 8855

dx.doi.org/10.1021/jp203220g | J. Phys. Chem. A 2011, 115, 8855–8859

The Journal of Physical Chemistry A

ARTICLE

malonic acid, and ferroin.10 Thus, we precisely know the total concentration of the catalyst C, but we do not know the concentration of HBrO3 nor the concentration of bromomalonic acid. To link these concentrations with the initial solution of reagents, we assume that the concentration of HBrO3 is proportional to the concentration of NaBrO3 used in experiment A ¼ ½HBrO3  ¼ η½NaBrO3  ¼ ηN

ð4Þ

and that the concentration of CHBr(COOH)2 is proportional to the product of [CH2(COOH)2] and [KBr] B ¼ ½CHBrðCOOHÞ2  ¼ ξ½CH2 ðCOOHÞ2 /½KBr ¼ ξM/K

ð5Þ

In eqs 4 and 5, K, M, and N denote concentrations of KBr, CH2(COOH)2, and NaBrO3, respectively. It is convenient to consider K, M, and N as parameters of the model because concentrations of these reagents are fixed when an experiment starts. Introducing the scaled concentration of the oxidized catalyst Figure 1. Time evolution of y within one cycle. The solid line represents the three-variable model and the dashed one the results of the twovariable model. Calculations were performed for the following conditions: [H2SO4] = 0.30 M, [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.35 M, [KBr] = 0.06 M, and C = 0.0017 M.

chemical oscillations do not seem to be affected by the presence of phospholipids. The surrounding lipid layer stabilizes a droplet against merging with the others. The experiments have demonstrated10 that if a droplet is not too small (its diameter is on the order of 1 mm or larger), then the period of homogeneous oscillations in such a droplet is the same as the period of oscillations in a test tube with a well-stirred solution with the same composition of reagents. Therefore, we can describe the time evolution of reagents in droplets with a simple model of the BZ reaction that does not take into account reactions in which lipids are involved. A simple, few-variable model can be derived from the Field
Kores
Noyes (FKN)14 reaction scheme. After some simplifications,15
17 the FKN reaction scheme can be reduced to the following kinetic equations ∂X ¼ k1 h0 AX
2k4 h0 X 2
k5 h0 XY þ k7 h0 AY ∂t

ð1Þ

∂Y k8 k9 B Z ¼q
k5 h0 XY
k7 h0 AY þ k13 B ∂t k
8 h0 ðC
ZÞ

ð2Þ

∂Z k8 k9 B Z ¼ 2k1 h0 AX
∂t k
8 h0 ðC
ZÞ

z ¼ Z=C and the scaled concentration of X defined as: x ¼ 2k1 ηX one transforms eq 3 into ∂z h0 N K/M z ¼ x
R ∂t C Ch0 ð1
zÞ

ð6Þ

where k8 k9 ξ k
8

R¼

If we introduce the scaled concentration of Y in the form y¼

Y k13 ξ

then eq 2 transforms into ∂y M/K z ¼ qβ
γh0 xy
γμh0 Ny þ M/K ∂t h0 1
z

ð7Þ

where

ð3Þ

The model variables X, Y, and Z denote concentrations of HBrO2, Br
, and Fe(phen)3+ 3 , respectively. In eqs 1
3, the symbols k+
i represent the rate constants of the corresponding reaction within the FKN reaction scheme.17 Here, h0 is the Hammett acidity function of the solution, and for the concentrations of sulfuric acid used in our experiments, it can be approximated as 1.3 times the concentration of H2SO4. The capital letters A, B, and C denote the concentrations of HBrO3, bromomalonic acid CHBr(COOH)2, and the total concentra3+ tion of the catalyst (C = [Fe(phen)2+ 3 ] + [Fe(phen)3 ]). We assume that the values A, B, and C remain constant during the observed evolution. In our experiments, BZ droplets were prepared with a water solution of sulfuric acid, sodium bromate, potassium bromide,

γ¼

k5 2k1 η

μ¼

k7 η 2k1 k7 η2 ¼ γ k5

β¼

k8 k9 k
8 k13

Finally, within such scaling, eq 1 can be written as ∂x RγE1 RγE1 μ ¼ E1 h0 Nx
E2 h0 x2
2 h0 Ny h0 xy þ 2 ∂t β β

ð8Þ

where E1 ¼ k1 η E2 ¼ 8856

k4 k4 ¼ E1 k1 η dx.doi.org/10.1021/jp203220g |J. Phys. Chem. A 2011, 115, 8855–8859

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Period as a function of concentration of different reagents. Circles show the experimental results, the solid line has been calculated using the three-variable model, and the dashed one plots results of the two-variable model. The periods were calculated for the following conditions: (a) [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.35 M, [KBr] = 0.06 M, and C = 0.0017 M; (b) [H2SO4] = 0.30 M, [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.35 M, and C = 0.0017 M; (c) [H2SO4] = 0.45 M, [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.35 M, and [KBr] = 0.06 M; (d) [H2SO4] = 0.30 M, [CH2(COOH)2] = 0.35 M, [KBr] = 0.06 M, and C = 0.0017 M.

Let us notice that all parameters of the model (R, β, γ, ɛ1, ɛ2, μ, and q) depend on the rate constants only, and they are not related to concentrations of any considered reagents. Therefore, the three-variable model based on eqs 6
8 can be directly used to describe the time evolution of a solution of reagents prepared with known concentrations of H2SO4, KBr, CH2(COOH)2, NaBrO3, and the catalyst. We think that it makes it superior to the Rovinsky
Zhabotinsky model of the BZ reaction where all model parameters have to be separately calculated for each experimental case. The system based on equations 6
8 can be made simpler if we assume that the relaxation of y is fast and the stationary value of y is immediately approached. Calculating the steady-state value of y from eq 7 and substituting it into eq 8, we obtain ∂x ¼ E1 h0 Nx
E2 h0 x2 ∂t   1 1 z x
μN þ q
2RE1 M/K β h0 1
z x þ μN

ð9Þ

The set of equations (eqs 6 and 9) with parameters R, β, ɛ1, ɛ2, μ, and q constitutes the two-variable model describing the BZ reaction in lipid-covered droplets.

3. EXPERIMENTAL RESULTS AND PARAMETER FITTING In the experiments, we measured the period of oscillations in a small droplet of BZ reagents surrounded by a phospholipid layer. The medium used was composed of two immiscible phases. One was the solution of a phospholipid L-R-phosphatidylcholine (Soy-20%, Avanti Polar Lipids, Inc.) in decane, prepared by dissolving 1 g of lipids in 50 mL of decane. The second phase was the aqueous solution of reagents of ferroin-catalyzed BZ reaction corresponding to an oscillating regime. The droplets were individually placed in small test tubes containing solutions of phospholipids.10 The experimental data presented below were obtained for droplets with diameters between 1 and 1.5 mm. We restricted the droplet size because for larger droplets, the probability that a nonhomogeneous structure spontaneously appears is large. The experiments were performed for over 60 different compositions of reagents. In the ferroin-catalyzed BZ 8857

dx.doi.org/10.1021/jp203220g |J. Phys. Chem. A 2011, 115, 8855–8859

The Journal of Physical Chemistry A

ARTICLE

model independently and performed optimization of its parameters for the same set of experimental results. For the twovariable model, the optimization procedure gives the best fit for R = 2.6  10
4, β = 200, ɛ1 = 4000, ɛ2 = 5800, μ = 2.1  10
5, and q = 0.88. The values of all parameters except R are quite close to those optimized for the three-variable model. The periods calculated with the two-variable model are shown by the dashed line in Figure 2. The agreement is not as good as that for the three-variable model, but still, most of the trends in period changes are correctly represented.

Figure 3. Nullclines for the two-variable model. The solid line corresponds to the x-variable (eq 9) and the dashed one to the z-variable (eq 6). The upper pair of curves represents a typical oscillatory system ([H2SO4] = 0.22 M, [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.35 M, [KBr] = 0.06 M, and C = 0.0017 M), and the lower pair describes an excitable medium ([H2SO4] = 0.02 M, [NaBrO3] = 0.45 M, [CH2(COOH)2] = 0.50 M, [KBr] = 0.168 M, and C = 0.0017 M).

reaction, changes in the concentration of reagents are indicated by changes in color; the solution is blue for a high concentration of oxidized catalyst and red when the catalyst is reduced. Therefore, information on the time evolution and period can be obtained from visual observation of the system. The time evolution of the BZ droplets was recorded with a digital video camera. Next, the movie was digitally processed, and the period of oscillations was calculated. The model parameters have been obtained from the experimental results for the period as a function of reagent concentrations. We used an optimization procedure that minimized differences between the observed period and the one calculated from model equations for a large set of experimental conditions. For the model based on eqs 6
8, the best agreement was obtained for R = 1.1  10
4, β = 200, γ = 6000, ɛ1 = 4300, ɛ2 = 8800, μ = 2.5  10
5, and q = 0.6. With these values of parameters and concentrations K, M, and N given in mol/L, the second is the time unit in eqs 6
8. The periods calculated from the three-variable model with such parameters are shown by the solid line in Figure 2 and compared with the experimental results (open cycles). It can be seen that this model gives a quite accurate description of the trends in period changes as the concentration of a selected reagent is varied. Using the parameters of the three-variable model, we can verify if the assumption on the stationarity of the y variable is correct. Figure 1 illustrates one particular case. The solid line shows the time evolution of y(t) predicted by the three-variable model. The dashed line plots the value of ys calculated from values of x(t) and z(t) within the assumption that ys corresponds to the steady state of eq 7. It can be seen that for almost all periods, the agreement is quite good and y(t) is roughly 10% larger than its estimation ys. Only around the excitation do the numbers differ by an order of magnitude. Having this difference in mind, one can expect that the parameters optimized for the three-variable model are not the best ones for the twovariable model based on eqs 6 and 9. We treated the two-variable

4. CONCLUSIONS In this paper, we have introduced two simple models of the BZ reaction in which concentrations of reagents used to initiate experiments appear as the input variables. The numerical values of model parameters have been fixed such that both models optimally predict the periods of oscillations observed in a large set of experiments performed at different conditions. We think that the models with the proposed values of parameters can be used in simple yet qualitatively correct simulations of the time evolution of BZ droplets surrounded by hydrocarbons. They can be also applied to predict experimental conditions at which the requested character of evolution appears. For example, the twovariable model exhibits excitability when C = 0.0017, the sodium bromate concentration is 0.45, the concentration of sulfuric acid is reduced to 0.02, and the product of concentrations of malonic acid and potasium bromide is increased to 0.084. The related changes in nullcline positions are illustrated in Figure 3. We think that the proposed models can form the basis for future studies on spatiotemporal effects in multiple droplet systems that lead to information processing operations. In order to describe interactions between different droplets, we have to generalize the models introducing the diffusion within a droplet and the transport of reagents in lipid layers and hydrocarbons. The variables y and z describe ionic reagents that do not diffuse in hydrocarbons. Moreover, the variable z describes a large complex, and its diffusion in water solution seems to be negligible if compared with the other reagents. Therefore, the generalization of the two-variable model seems quite simple because it is sufficient to estimate the diffusion constant of HBrO2 in the medium and its migration through the lipid layer on the basis of experimental results. Within the threevariable model, one should additionally consider diffusion of Br
in the water phase. We will return to the problem of estimation of mobility parameters in our subsequent paper. ’ ACKNOWLEDGMENT The research was supported by the NEUNEU project sponsored by the European Community within FP7-ICT-2009-4 ICT-4-8.3 - FET Proactive 3: Biochemistry-based Information Technology (CHEM-IT) program. ’ REFERENCES (1) Hjelmfelt, A.; Weinbergret, E. D.; Ross, J. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 10983–10987. (2) Buisman, H. J.; ten Eikelder, H. M. M.; Hilbers, P. A. J.; Liekens, A. M. L. Artificial Life 2009, 15, 1–15. (3) Adamatzky, A.; De Lacy Costello,B.; Asai, T. Reaction
Diffusion Computers: Elsevier Science: Amsterdam, The Netherlands, 2005. (4) Gorecki, J.; Gorecka, J. N. Computing in Geometrical Constrained Excitable Chemical Systems. In Encyclopedia of Complexity and Systems Science, Meyers, R. A. ed., Springer-Verlag: New York, 2009. 8858

dx.doi.org/10.1021/jp203220g |J. Phys. Chem. A 2011, 115, 8855–8859

The Journal of Physical Chemistry A

ARTICLE

(5) Tero, A.; Kobayashi, R.; Nakagaki, T. Physica A 2006, 363, 115– 119. (6) Kuhnert, L.; Agladze, K. I.; Krinsky, V. I. Nature 1989, 337, 244– 247. (7) Thakore, V.; Behal, A.; Molnar, P.; Leistritz, D. C.; Hickman, J. J. J. Comput. Theor. Nanosci. 2008, 5, 2164–2169. (8) Vanag, V. K.; Epstein, I. R. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 14635–14638. (9) Toiya, M.; Vanag, V. K.; Epstein, I. R. Angew. Chem., Int. Ed. 2008, 47, 7753–7755. (10) Szymanski, J; Gorecka, J. N.; Igarashi, Y.; Gizynski, K.; Gorecki, J.; Zauner, K. P.; de Planque, M Int. J. Unconv. Comput. 2011in press. (11) Adamatzky, A.; De Lacy Costello, B.; Bull, L.; Holley, J. Isr. J. Chem. 2011, 51, 56–66. (12) Holley, J.; Adamatzky, A.; Bull, L.; De Lacy Costello, B.; Jahan, I. Nano Commun. Networks 201110.1016/j.nancom.2011.02.002. (13) Brandstadter, H.; Braune, M.; Schebesch, I.; Engel, H. Chem. Phys. Lett. 2000, 323, 145–154. (14) Field, R. J.; Koros, E.; Noyes, R. M. J. Am. Chem. Soc. 1972, 94, 8649–8664. (15) Rovinsky, A. B.; Zhabotinsky, A. M. J. Phys. Chem. 1984, 88, 6081–6084. (16) Rovinsky, A. B. J. Phys. Chem. 1986, 90, 217–219. (17) Sielewiesiuk, J.; Gorecki, J. J. Phys. Chem. A 2001, 105, 8189– 8195.

8859

dx.doi.org/10.1021/jp203220g |J. Phys. Chem. A 2011, 115, 8855–8859