Spiral, Core-defect and Wave break in a modified Oregonator Model

Parvej Khan and Sumana Dutta ¹¹1sumana@iitg.ac.in Department of Chemistry, Indian Institute of Technology Guwahati,
Guwahati 781039, INDIA.

Abstract

Target waves and spiral waves were discovered in the Belousov-Zhabotinsky (BZ) reaction around 50 years ago. Many biological systems demonstrate such rotating spiral patterns. Spiral waves are widely encountered in the glycolytic activity of yeast and cardiac and neuronal tissues. In the cardiac system, the spiral waves and their three dimensional counterparts, the scroll waves are the cause of arrhythmia. These waves are also widely studied in the experimental BZ reaction as a table-top model of the cardiac system. To understand the underlying physics or chemistry of the system, the Oregonator model is considered in correlation to the BZ reaction system. In this study, starting from the same FKN mechanism, we slightly modify the original Oregonator model. We present a three-variable model that can be called a modified Oregonator model in which we observe the presence of the spiral, spiral defect, and spiral breakup by tuning the parameters.

Belousov-Zhabotinsky reaction

|

Spiral Breakup

|

Modified Oregonator model

|

…

Patterns like spirals and scrolls are abundant in nature. A large number of systems show the formation of rotating spiral waves win . From cardiac tissue to BZ reaction-in a wide range of biological, physiological, physical, and chemical systems can be noticed. One of the main ambitions or goals is to understand the electrical activity of the human cardiac system through the study of a spiral and its three-dimensional counterpart jalife .

In most of the observed systems, the spiral rotates outwardly from its core, which is known as a pacemaker copt ame glyco . In connection to the cardiac system, a pair of spiral waves in a BZ system are investigated in many ways. The interaction of multiple counter-rotating spiral pairs are being studied by Kalita et al. both experimentally as well as in numerical simulation hk_int . Spiral and scroll dynamics are being tried to control by different kinds of external gradients like thermal gradients, electric field gradients etc. luther ,mix . Spiral wave properties are also being investigated with variation of excitability factor Dhriti_spiral .

One of the very interesting as well as important phenomena is spiral break up or wave break up. Spiral breakup is being studied extensively in a computer simulation using CGLE and FHN models turb_bu , though we have limited examples in experiments. Sustainable spiral waves in the human heart are comparable with cardiac arrhythmia and the wave breakup is related to more dangerous situations like ventricular fibrillation bu_vfibri . So the dynamics of the spiral wave and mechanism of wave breakup is very necessary to study both in experiments also in theory.

Belousov-Zhabotinsky (BZ) reaction is among the most studied oscillatory chemical reactions. It can generate sustainable chemical waves like spirals and scrolls winfree . Though there are differences in heterogeneity, conduction velocity, etc., these waves are excitable as the electrical signal of our cardiac system. Theoretical treatment of the mechanistically complex BZ reaction is done with the FKN mechanism and a well-established Oregonator model is often used to generate spiral and scroll waves to study.

In this manuscript, we describe a newly developed model derived from the FKN mechanism, which may also be called as modified Oregonator model. We investigated the model through numerical simulations. We report the formation of spiral, annihilation of spirals and also the exciting wave break phenomena upon parameter tuning. Previously, almost in all studies, specific models are generally used to describe some specific phenomena but here we show almost all kinds of spiral dynamics with this single four-variable reaction-diffusion model.

I Model

I.1 Mechanism for BZ reaction

The important steps in the mechanism of the BZ reaction proposed by Field-Körös-Noyes(FKN) spiral_scott , used for development of most numerical models, is given as $-$

$\displaystyle\text{BrO}_{3}^{-}+\text{Br}^{-}+2\text{H}^{+}$	$\displaystyle\xrightarrow{k_{1}}$	$\displaystyle\text{HBrO}_{2}+\text{HOBr}$	(1)
$\displaystyle\text{HBrO}_{2}+\text{Br}^{-}+\text{H}^{+}$	$\displaystyle\xrightarrow{k_{2}}$	$\displaystyle 2\text{HOBr}$	(2)
$\displaystyle\text{BrO}_{3}^{-}+\text{HBrO}_{2}+\text{H}^{+}+2\text{M}^{2+}$	$\displaystyle\xrightarrow{k_{3}}$	$\displaystyle 2\text{HBrO}_{2}+2\text{M}^{3+}+\text{H}_{2}\text{O}$	(3)
$\displaystyle 2\text{HBrO}_{2}$	$\displaystyle\xrightarrow{k_{4}}$	$\displaystyle\text{BrO}_{3}^{-}+\text{HOBr}+\text{H}^{+}$	(4)
$\displaystyle 2\text{M}^{3+}+\text{MA}+\text{BrMA}$	$\displaystyle\xrightarrow{k_{5}}$	$\displaystyle f\text{Br}^{-}+2\text{M}^{2+}+\text{Other products}.$	(5)

The FKN mechanism involves three main processes: the autocatalysis of HBrO₂ (step 3); inhibition by bromide ion (steps 1 and 2), and resetting the cycle (4 and 5). The values of the rate constants are given as $k_{1}$ = 2 M^-3s^-1, $k_{2}$ = 1.25 $\times$ 10⁵ M^-2s^-1, $k_{3}$ = 45 M^-2s^-1, $k_{4}$ = 1.4 $\times$ 10³ M^-1s^-1, $k_{5}$ = 6.5 $\times$ 10^-3 M^-1s^-1 spiral_scott .

Expressing the reactions in a more general form gives-

$\displaystyle\textbf{A}+\textbf{W}+2\textbf{H}$	$\displaystyle\xrightarrow{k_{1}}$	$\displaystyle\textbf{U}+\textbf{P}$	(6)
$\displaystyle\textbf{U}+\textbf{W}+\textbf{H}$	$\displaystyle\xrightarrow{k_{2}}$	$\displaystyle 2\textbf{P}$	(7)
$\displaystyle\textbf{A}+\textbf{U}+\textbf{H}$	$\displaystyle\xrightarrow{k_{3}}$	$\displaystyle 2\textbf{U}+2\textbf{V}$	(8)
$\displaystyle 2\textbf{U}$	$\displaystyle\xrightarrow{k_{4}}$	$\displaystyle\textbf{A}+\textbf{P}+\textbf{H}$	(9)
$\displaystyle\textbf{V}+\textbf{B}$	$\displaystyle\xrightarrow{k_{5}}$	$\displaystyle{\frac{1}{2}}f\textbf{W}$	(10)

This system of equations constitute the Modified Oregonator model, where comparison with the FKN would yield A = BrO ${}_{3}^{-}$ , W= Br^-, H= H⁺, U= HBrO₂, P = HOBr, V = oxidized form of catalyst, and B = MA+BrMA (malonic acid and its derivatives). This scheme is able to give, by simple linear combination of equations, the overall reaction, $f\textbf{A}+4\textbf{B}+(3f+1)\textbf{H}\rightarrow(3f+1)\textbf{P}$ , with no net production or destruction of U, V and W. Hence, they are considered as the intermediates, whose dynamics is of interest to us.

A and B are major reactants, and their concentrations are treated as constants, a and b. As the reaction is carried out in an acidic environment, we have an excess of hydrogen ions at any time, and H is also considered as a constant parameter.

The kinetic behavior of the intermediates can now be written from the above model as-

$\displaystyle\frac{dU}{d\tau}$	$\displaystyle=$	$\displaystyle k_{1}aWH^{2}-k_{2}UWH+k_{3}bUH-2k_{4}U^{2}\$	(11)
$\displaystyle\frac{dV}{d\tau}$	$\displaystyle=$	$\displaystyle 2k_{3}aUH-k_{5}bV\$	(12)
$\displaystyle\frac{dW}{d\tau}$	$\displaystyle=$	$\displaystyle-k_{1}aWH^{2}-k_{2}UWH+\frac{1}{2}k_{5}fbV\$	(13)

Here $U$ , $V$ , $W$ , and $H$ are the concentrations of the respective chemical species in units of mol L^-1 and $\tau$ is time in units of s.

I.2 Dimensionless form

We further non-dimensionalize the equations in the above sequence (11-13), using the following scaling.

\displaystyle u=\frac{k_{3}}{k_{5}}\frac{a}{b}U,\hskip 14.22636ptv=\frac{k_{1}k_{3}}{k_{2}k_{5}}\frac{a^{2}}{b}V,\hskip 14.22636ptw=\frac{k_{2}}{k_{3}}\frac{1}{a}W,\hskip 14.22636ptt=k_{5}b\tau,

where $u$ , $v$ , and $w$ , are the dimensionless variables, giving us the set of non-dimensional kinetic equations

$\displaystyle\epsilon\frac{du}{dt}$	$\displaystyle=$	$\displaystyle\frac{1}{2}wh^{2}-uwh+uh-qu^{2}\$	(14)
$\displaystyle\frac{dv}{dt}$	$\displaystyle=$	$\displaystyle uh-v$	(15)
$\displaystyle\epsilon^{\mathrm{{}^{\prime}}}\frac{dw}{dt}$	$\displaystyle=$	$\displaystyle-qwh^{2}-2quwh+2qfv\$	(16)

\noindent\text{where }h=\frac{2k_{1}k_{3}}{k_{2}k_{5}}\frac{a^{2}}{b}H=\frac{1}{h_{0}}H,\text{is the non-dimensional parametrized form of the }

hydrogen ion concentration, H. In the above set of equations, the parameters are given by,

\displaystyle\epsilon=\frac{1}{h_{0}}\frac{k_{5}}{k_{3}}\frac{b}{a},\hskip 14.22636pt\epsilon^{\mathrm{{}^{\prime}}}=\frac{1}{h_{0}^{2}}\frac{2k_{4}k_{5}}{k_{2}k_{3}}\frac{b}{a},\hskip 14.22636pt\text{and }\hskip 8.5359ptq=\frac{2k_{1}k_{4}}{k_{2}k_{3}}.\

The values of $a$ and $b$ are considered to meet usual experimental conditions oreg , and can be varied to give different values of the parameters.

Now we have a modified version of the three-variable Oregonator model as equations (14-16).

II Simplified two variable model

The BZ-system is traditionally studied by the Oregonator model which is derived from the FKN (Field- Körös -Noyes) mechanism FKN1972 of the reaction. The Oregonator model is a three-variable system, that is often simplified to a two-variable activator inhibitor model, where bromous acid (HBrO₂) is the activator, $u$ , and the oxidized form of the catalyst,(Fe³⁺) is the inhibitor, $v$ oreg .

	$\displaystyle\frac{du}{dt}$	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon}\left(u(1-u)-\frac{fv(u-q)}{(q+u)}\right)$		(17)
	$\displaystyle\frac{dv}{dt}$	$\displaystyle=$	$\displaystyle u-v$		(18)

However, the hydrogen ion concentration does not appear in either the two- or three- variable models of the Oregonator, as it gets embedded in the kinetic parameters $f$ , $q$ and $\epsilon$ . Our modified model described above considers hydrogen ion as an important parameter. The three variable model described above can be more simplified.

As, $\epsilon^{\mathrm{{}^{\prime}}}$ $<<$ $\epsilon$ , steady state approximation can be employed on $w$ oreg , giving us

w=\frac{2fv}{h(h+2u)}

Substituting this value of $w$ in the above equations [14-16], we have:

	$\displaystyle\epsilon\frac{du}{dt}$	$\displaystyle=$	$\displaystyle fv\frac{h-2u}{h+2u}+u(h-qu)$		(19)
	$\displaystyle\frac{dv}{dt}$	$\displaystyle=$	$\displaystyle uh-v\$		(20)

In the presence of diffusion, the modified Oregonator model takes the following form.

	$\displaystyle\epsilon\frac{du}{dt}$	$\displaystyle=$	$\displaystyle fv\frac{h-2u}{h+2u}+u(h-qu)+D_{u}\nabla^{2}u$		(21)
	$\displaystyle\frac{dv}{dt}$	$\displaystyle=$	$\displaystyle uh-v+D_{v}\nabla^{2}v$		(22)

where $D_{u}$ and $D_{v}$ are the diffusion coefficients of $u$ and $v$ respectively.

Refer to caption — Figure 1: Stable spiral. (a)-(c) Snapshots at 324.2 t.u.,327.1 t.u., and 330.3 t.u., respectively. The star marks on the spiral wave-arm show the outwardly movement of the wave. (d) Time variation of the three variables $u$ as solid black curve, $v$ the dashed red curve and $w$ is the dash dotted blue curve (same curve types are used across the images) (e) Time space plot between 150 to 200 t.u. along yellow arrow shown in (a). Color key for the snapshots and timespace plots is depicted between (a) and (b), and is similar across the figures.

III Generation of spiral waves

In order to study spiral wave activity in a spatially extended system, the BZ reaction is often carried out in an unstirred gel medium where diffusion plays a major role in modifying the dynamics of the system. The gel however, rules out any convection in the medium. For such a system, our three variable reaction-diffusion equations are given as:

$\displaystyle\frac{du}{dt}$	$\displaystyle=$	$\displaystyle\frac{1}{\epsilon}\left[\frac{1}{2}wh^{2}-uwh+uh-qu^{2}\right]+D_{u}\nabla^{2}u\$	(23)
$\displaystyle\frac{dv}{dt}$	$\displaystyle=$	$\displaystyle\left[uh-v\right]+D_{v}\nabla^{2}v\$	(24)
$\displaystyle\frac{dw}{dt}$	$\displaystyle=$	$\displaystyle\frac{1}{\epsilon^{\mathrm{{}^{\prime}}}}\left[-qwh^{2}-2quwh+2qfv\right]+D_{w}\nabla^{2}w\$	(25)

where, $D_{u}$ , $D_{v}$ , $D_{w}$ are the diffusion coefficients of $u$ , $v$ and $w$ , respectively.

In order to numerically integrate the equations (23–25), zero flux boundary conditions are applied across all four boundaries. The system parameters are chosen as $\epsilon=0.035$ , $\epsilon^{\mathrm{{}^{\prime}}}=0.008$ , $f=0.9$ , $q=0.01$ . The time-step and space-step were chosen as $\Delta t=0.015$ t.u. and $\Delta x=0.45$ s.u., respectively, which translates to 0.814 ms, and 4.06 $\mu$ m. These values of $\Delta t$ and $\Delta x$ are convergent and yield best resolution of patterns in our parameter range, and also can be varied by around 20 $\%$ , without causing any major changes to the nature of the waves, or the time period and wavelength of the corresponding spirals. Simulations were carried out in a $300\times 300$ (135 s.u. $\times$ 135 s.u.) grid. Euler’s method is used for the integration, and a central difference method helped to calculate the diffusion term. Stable rotating spirals, with non-meandering circular cores are observed in a range of parameters. For our analysis, we vary the parameters from a base condition, at which the system supports an outwardly rotating spiral. The diffusion constants of $D_{u}=0.3$ , $D_{v}=0.17$ , $D_{w}=0.34$ , with $h=$ 0.3. The values of the diffusion coefficients are taken in a ratio so as to match their actual expected values depending on their molecular weights (molecular weight of $u$ or HBrO₂ is 113 g mol^-1, of $v$ or ferroin is 596 g mol^-1 and that of $w$ or bromide is 80 g mol^-1.). The initial values of $u$ , $v$ and $w$ are taken to be zero across the entire space, as they are intermediates of the BZ reaction, except for a thin area signifying the initial wave-front. A plane wave is initiated in the middle of the square area, by taking three stripes of width $2\Delta x$ and length $60\Delta x$ having non-zero values of the variables, emulating the front and back of the wave. The straight plane wave initiated, curled up to form a pair of counter-rotating spirals 1.

III.1 Effect of hydrogen ion concentration

Earlier studies on the BZ system and the Oregonator model had shown that changing hydrogen ion concentration had a major effect on the excitability of the system, and hence modified the wave nature Dhriti_spiral . However, in those studies, the effect of the H ${}^{+}-$ ion concentration could only be considered implicitly as a modification of the system parameter, $\epsilon$ , or an added acidity factor. Our modified Oregonator model enables us to explicitly vary the initial hydrogen ion concentration to predict the effect of the same on the nature of the system.

In order to explore the effect of $h$ , we start with the base condition, and vary the value of $h$ . Increasing the value of $h$ upto 0.55, stable spirals are obtained, beyond which the system has no solution. Decreasing the value of $h$ however yields interesting results.

III.1.1 Drifting Spirals

As $h$ is decreased below 0.29, the spiral starts to get unstable. Now, the spiral tips become anisotropic, and starts drifting away from their initial positions.

III.1.2 Spiral breakup

At still lower $h$ values, the entire system starts getting unstable and wave-break starts initiating far from the core, finally leading to the formation of multiple spirals in the system. The current scenario is reminiscent of the far-field breakup that has been observed in some earlier studies markusbar .

III.1.3 Target pattern

At $h\leq$ 0.17, there is no more formation of spiral. The straight plane wave that we initiate does curl up at its ends, but before they turn further to form a pair of counter-rotating spirals, the two ends touch and form target waves.

III.1.4 Chaotic Oscillations and Oscillation Death

Before we reach oscillation death, at values of $h<0.1$ , it passes through a fleeting region of chaotic oscillations, that is distinct from the spiral breakup regions.

For the base condition, a variation in $h$ , generates a 1-D phase diagram [2], showing all the above phenomenon. This suggests that merely changing the hydrogen ion concentration can generate a plethora of phenomena.

Our two variable model described above (Eqs. 21, 22) are also capable of generating spiral waves numerically. The presence of $h$ parameter helps to study the effect of hydrogen ions on the dynamics and wave properties of spirals in a BZ reaction-diffusion system.

IV Discussion and Conclusion

We report the formation of spiral waves, core defect, annihilation of spiral pairs, and the breakup of spirals. Spiral waves have both chemical and biological significance. These are abundant in cardiac tissues, catalytic surfaces (CO oxidation on Pt surface), and aggregation of amoebae. Spiral breakup near the core and outside the core is explained in the paper of Bar et. al. They took a more general Barkley model considering a very fast diffusion of the activator and delay in the production of inhibitor. We obtained both the spiral pair break-up scenario, where – (1) break up starts at the core and then spreads over space, and (2) break up starts away from the core leaving the core greatly unaffected. Many other papers concerned with the breakup dynamics of spirals also considered the FHN model or, the more general Barkley model with a faster diffusion of the activator 11 ; 13 ; 14 ; 15 . This current work is different in the sense that it uses a new model and gives a parametric regime of different kinds of spiral dynamics as well as the breakup phenomena. We also considered the diffusion coefficients of inhibitors like Br^- and Fe³⁺ are instead of neglect them like it has been done in other spiral breakup studies with Barkley or FHN model.

References

(1) A. T. Winfree, Spiral Waves of Chemical Activity,Science 175, 634 (1972).
(2) J. M. Davidenko, A. V. Pertsov, R. Salomonsz, W. Baxter,and J. Jalife, Nature 355,349,(1992).
(3) S. Jakubith, H. Rotermund, W. Engel, A. Von Oertzen, and G. Ertl, Phys. Rev. Let. 65, 24, 3013(1990).
(4) F. Siegert and J. C. Weijer,Physica D: Nonlinear Phenomena 49,224(1991).
(5) T. Mair and S. C. Müller, Journal of Biological Chemistry 271, 627 (1996).
(6) H. Kalita and S. Dutta, Phys. Rev. E 105, 054213 (2022)
(7) D. Mahanta, N.P. Das and S. Dutta, Phys. Rev. E 97, 022206 (2018)
(8) Luther et al. , Nature, 475, 235-239 (2011).
(9) Das, Nirmali Prabha and Dutta, Sumana, Physical Review E 96, 022206 (2017).
(10) M. B¨ar and M. Eiswirth, Phys. Rev. E 48, R1635 (1993).
(11) A. V. Panfilov, Chaos 8, 57 (1998).
(12) Winfree, AT, RoyalSociety of Chemistry 9, 38-46(1974).
(13) S. K. Scott, Oscillations, Waves, and Chaos in Chemical Kinetics, Oxford University Press, 1994.
(14) R. M. Noyes, R. Field, and E. Körös, J. Am. Chem. Soc. 94 (1972) 1394-1395.
(15) R. J. Field and R. M. Noyes, J. Chem. Phys. 60 (1974) 1877-84.
(16) M. Bär, L. Brusch, New. J. Phys. 6 (2004) 5.
(17) J. Taboada, A. P. Mu˜nuzuri, V. P´erez-Mu˜nuzuri, M. G´omez-Gesteira and V. P´erez-Villar, Chaos 4,519 (1994).
(18) L. Q. Zhou and Q. Ouyang, Phys. Rev. Let. 85, 1650 (2000).
(19) M. B¨ar and M. Eiswirth, Phys. Rev. E 48, R1635 (1993).
(20) J. Yang, F. Xie, Z. Qu and A. Garfinkel, Phys. Rev. Let 91, 148302 (2003).