Mathematical modelling of the vitamin C clock reaction

Chemical clock reactions are characterized by a relatively long induction period followed by a rapid ‘switchover’ during which the concentration of a clock chemical rises rapidly. In addition to their interest in chemistry education, these reactions are relevant to industrial and biochemical applications. A substrate-depletive, non-autocatalytic clock reaction involving household chemicals (vitamin C, iodine, hydrogen peroxide and starch) is modelled mathematically via a system of nonlinear ordinary differential equations. Following dimensional analysis, the model is analysed in the phase plane and via matched asymptotic expansions. Asymptotic approximations are found to agree closely with numerical solutions in the appropriate time regions. Asymptotic analysis also yields an approximate formula for the dependence of switchover time on initial concentrations and the rate of the slow reaction. This formula is tested via ‘kitchen sink chemistry’ experiments, and is found to enable a good fit to experimental series varying in initial concentrations of both iodine and vitamin C. The vitamin C clock reaction provides an accessible model system for mathematical chemistry.


Introduction
'Clock reactions' encompass many different chemical processes in which, following mixing of the reactants, a long induction period of repeatable duration occurs, followed by a rapid visible change. These reactions have been studied for over a century, with early examples including the work of H. Landolt on the sulphiteiodate reaction in the 1880s-for review, see Horváth & Nagypál [1], and the work of G. Forbes et al. [2] on G. Vortmann's thiosulphate -arsenite/arsenate reaction [3]. The origin of the term 'clock' is unclear, however Richards and Loomis in 1927 [4] referred to '. . .the long familiar iodine "clock" depending upon the reduction of potassium iodate by sulfurous acid. . .', suggesting that the term had already been in common use for some time.
Conway's 1940 article [5] referred to both reactions as common tools in the classroom for teaching the Law of Mass Action. An early example of detailed mathematical modelling of multi-step reactions via systems of nonlinear differential equations was given by Chien [6]; this work exhibited a number of analytical techniques including solving a Riccatitype equation resulting from quadratic reaction kinetics (a step we will find useful in what follows). Anderson [7] demonstrated the computer solution of chemical kinetics problems, including the iodine clock reaction. Further review and history of the field can be found in [1,8].
Through dynamical systems analysis and the method of matched asymptotic expansions, Billingham & Needham [8] identified and studied two clock reaction mechanisms, both alone and in combination. The first mechanism is autocatalysis: the reaction producing the clock chemical is catalysed by the clock chemical itself, for example, the reaction P þ 2B ! 3B, where B is the clock chemical and P is a precursor. This type of reaction is characterized by nonlinear kinetics; in the above case, the reaction rate for production of B would be proportional to pb 2 , where lower case letters correspond to chemical concentrations.
The second mechanism identified was inhibition. An inhibitor chemical C removes the clock chemical, e.g. B þ C ! D, keeping the concentration of B low. Simultaneously, clock chemical is produced upstream from a precursor, e.g. via the reaction P ! B. Provided that the supply of precursor is sufficiently large relative to the initial concentration of inhibitor, the inhibitor will eventually be depleted, allowing the clock chemical concentration to rise. Billingham & Needham [8] then applied this framework to the iodate-arsenous acid system, characterizing the process as involving a combination of inhibition (B þ C ! 2A) of the clock chemical and indirect autocatalysis in one of the reactants A. The catalysis is indirect because it is through the combination of this reaction with (P þ 5A ! 3B) results in A effectively catalysing its own production (5A produce 6A). For the iodate-arsenous acid system, the clock chemical B is iodine (I 2 ), the inhibitor C is arsenite ion AsO 3À 3 , the precursor P is iodate (IO À 3 ) and the other reactant A is the iodide ion (I 2 ). Further mathematical development can be found in [9,10], application of this approach to systems of industrial importance [11] and biochemistry [12]. Further discussion on the appropriateness of the term 'clock reaction' for various processes involving induction periods can be found in [13].
The system we will consider was described by Wright [14,15] as a variation on the iodine clock reaction requiring only safe household chemicals. A combination of iodine (B) and iodide ions (A) are supplied initially. In the presence of hydrogen peroxide, iodide ions are converted to iodine molecules, (1:1) The rate of this reaction has been considered by a number of authors [16][17][18], who have observed that the rate is linear in both I 2 and H 2 O 2 concentrations due to a rate-limiting step involving nucleophilic attack of I 2 on H 2 O 2 . However, these studies have worked with hydrogen peroxide concentrations which are similar to the other substrates, whereas we will work with hydrogen peroxide in great excess. In consequence, the H 2 O 2 concentration may be omitted from the model, and a one-step reaction with rate that is quadratic in I 2 concentration (as would be expected from applying the law of mass action) is found to match well to experimental results. We therefore have in simplified form, Starch in solution appears blue in the presence of iodine. There is no separate precursor P in this reaction. The inhibitor C is ascorbic acid-i.e. vitamin C-which converts iodine to iodide (alongside producing equal quantities of H þ ions) via the reaction, with initial conditions based on the initial concentrations of A, B and C, respectively.
Inspecting equations (2.1) and (2.2), it is clear that the total concentration of iodine atoms a(t) þ 2b(t) is conserved. Denoting the initial concentration by m 0 ¼ a 0 þ 2b 0 , we then have, which leads to the two-variable system, and dc dt ¼ Àk 1 bc: Non-dimensionalizing the system with the scalings, leads to the dimensionless initial-value problem, db dt ¼ Àbg þ er(1 À 2b) 2 (2:9) and dg dt ¼ Àrbg, where the dimensionless parameters r ¼ m 0 /c 0 and e ¼ k 0 /k 1 . A key feature of the dynamics is that because the reaction producing clock chemical is much slower than the inhibitory reaction, the latter ratio is very small, i.e. e ( 1. The ratio r will be assumed to be order 1; it will also turn out to be that rf :¼ b 0 /c 0 , 1 in order for the vitamin C supply to survive the initial transient.

Qualitative analysis of the dynamics
Significant insight into the induction period of the system (2.9) and (2.10) can be obtained through an informal quasi-steady analysis, that is to exploit the fact that the clock chemical B is slowly varying during this period. If db/dt % 0 in equation (2.9), then it follows that bg % er(1 2 2b) 2 . Substituting this expression into equation (2.10) for inhibitor depletion leads to, dg dt % Àer 2 (1 À 2b) 2 : Given that clock chemical concentration is small during the induction period, b % 0 and so equation (3.1) can be simplified as dg/dt % 2er 2 , therefore g % c 1 2 er 2 t, where c 1 is a constant which we will determine in §4.2. This expression shows that the length of the induction period scales with (er 2 ) 21 in dimensionless variables.
The system (2.9) and (2.10) has the unique equilibrium (b, g) ¼ (1/2, 0) which represents the longterm fate of the system (zero inhibitor, all reactants converted to clock chemical). The eigenvalues of royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181367 the linearized system at this point are 0, with eigenvector (1, 0) and 2r/2, with eigenvector (1, r). The zero eigenvalue indicates the slow manifold {(s, 0) : s [ R} which we will find the system approaches after the induction period is complete. The negative eigenvalue indicates a stable manifold {(s, rs) : s [ R} which lies outside the physically reasonable region of interest 0 b 1/2, 0 g.

Asymptotic analysis
The smallness of the reaction rate ratio e enables an approximate solution to be sought via matched asymptotic expansions, in a similar manner to Billingham & Needham [8,9]. A visual summary of the asymptotic regions in the phase plane is shown in figure 1. We will assume throughout that the initial I 2 and vitamin C concentrations are such that rf , 1.

Initial adjustment (region I)
The first region that can be identified is where the independent variable t and both dependent variables b, g are order 1. Seeking a solution of the form, we find at leading order, is therefore constant, leading to the separable ordinary differential equation, It follows that As t ! 1 we observe that (b 0 (t), g 0 (t)) ! (0, 1 2 rf). This behaviour is the initial transient through which the system rapidly adjusts to quasi-equilibrium. During this interval, approximately rf ¼ b 0 /c 0 of the initial quantity of inhibitor is consumed. The next order terms in the dynamics may be found through straightforward but lengthy manipulations which do not reveal any significant further insight.

Induction (region II)
, the system takes the form, and dĝ dT ¼ Àrbĝ: Again seeking asymptotic expansionsb ¼b 0 þ eb 1 þ Á Á Á andĝ ¼ĝ 0 þ eĝ 1 þ Á Á Á we have the leading order problem, which has solutions of the form, Taking T ! 0 and matching to the region I solution as t ! 1 shows that the constant c 1 ¼ 1 2 rf. In the original dimensionless variables, the leading order solution in region II is therefore, The solution (4.12) becomes non-uniform as t ! (1 2 rf)/(r 2 e), which corresponds to the end of the induction period. In dimensional variables, the 'switchover' time is therefore characterized by, (4:13) We will return to equation (4.13) to interpret the experimental results in §5.

Long-term state (region IV)
We now turn our attention to the asymptotic dynamics in the long-term state of the system after the switchover, in which b(t) is again O(1) and t ¼ O(e 21 ). The corner region III between II and IV will be addressed later. Denoting b(T) ¼ b(e À1 T) and g(T) ¼ g(e À1 T), the system takes the form, and e d g dT ¼ Àr b g: (4:15) Again substituting expansions of the form b ¼ b 0 þ e b 1 þ Á Á Á and g ¼ g 0 þ e g 1 þ Á Á Á yields at leading order, b 0 g 0 ¼ 0, (4:16) hence g 0 ¼ 0, and moreover at O(e n ), royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181367 Hence g 0 ¼ Á Á Á ¼ g nÀ1 ¼ 0 implies that g n ¼ 0. By induction it follows that g n ¼ 0 for all n. Therefore, once b ¼ O(1) and T ¼ O(e 21 ), then g is zero at all orders in e. Taking g % 0 then yields, (note that we are now working with b rather than just the leading order term) which has solution, b(T) t ¼ e m (et À [r À2 À r À1 f]) and rescaled variables b ¼ e n b and g ¼ e p g we find that the most structured balance occurs when m ¼ n ¼ p ¼ 21/2.
The rescaled system is then, þ Á Á Á and similarly for g, the leading order problem is, The system (4.24) can be rearranged to yield the Riccati equation [8], where c 3 is a constant. Seeking a solution of the form g 0 ¼ u 0 =u yields the separable equation u 00 ¼ À(r 2 t þ c 3 )u 0 , and hence, The solution to the system can then be expressed in terms of the error function as and (4:28) The unknown constant c 4 can be fixed by considering the asymptotic form of g 0 as t ! À1. Using the asymptotic form erf(x) 1 À e Àx 2 x À1 p À1=2 (1 þ O(x À2 )) and for brevity defining f( t) :¼ ( )) as t ! À1: (4:29) Hence, Since f( t) ! À1 as t ! À1 we deduce that for g 0 ( t) to tend to a non-zero limit we must have that c 4 ¼ 1. In this case, the asymptotic behaviour is then ) as t ! À1: (4:31) In the region II variables, we then have (4:32) Matching to the region II solution then yields c 3 ¼ 0.
The approximate solutions in region III in the original variables are therefore, (4:34)

Comparison of asymptotic and numerical approximations
A comparison of the asymptotic approximations against a numerical solution of the system (calculated with the stiff solver ode15s, Matlab, Mathworks) is shown in figure 4, with the reaction rate ratio e taken as 0.001, reactants ratio r ¼ 2 and iodide : iodine ratio f ¼ 0.2 (these values are chosen arbitrarily). The region I solution follows the numerical approximation very closely up to around t ¼ 5. The region II solution then follows the numerical solution to within a few per cent up to around t ¼ 100. The region III solution then follows the numerical solution closely around the dimensionless switchover time e 21 [r 22 2 r 21 f ] ¼ 150, up to around t ¼ 200. Finally, the agreement between the region IV solution for t . 200 is excellent, as would be expected from the smaller-than-algebraic error in equations (4.41) and (4.42).

Experiments
In this section, we present the results of some 'kitchen sink' experiments, conducted with readily available chemicals. The protocol is essentially as described by Wright [14] (Procedure D. Reaction Using Kitchen Measuring Ware) only with 3% Lugol's iodine instead of tincture of iodine 2%, and varying the quantities of water and iodine used.
royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181367 8 Vitamin C stock solution was prepared by mixing a 1000 mg vitamin C tablet with 30-120 ml water. Solution 'A' is prepared from 5 ml (1 teaspoon) of vitamin C stock solution, between 2.5 and 10 ml of 3% w/v Lugol's iodine, and 60 ml (4 tablespoons) of warm tap water at 408C. Solution 'B' is prepared from 15 ml of 3% hydrogen peroxide, 1 level teaspoon of powdered laundry starch and 60 ml of warm water. Solutions A and B are added in a beaker, at which point a timer is started, and the mixture is stirred manually for 5 s. The timer is stopped when the colour change from clear to blue occurs.
Lugol's iodine is formulated as a 1 : 2 mixture of iodine (I 2 ) and potassium iodide (KI) [20]. Based on molar masses of 126.9 g mol 21 for iodine and 166 g mol 21 for potassium iodide, Lugol's 3% w/v solution contains 1 þ 2 Â 126.9/166 ¼ 2.5289 g/100 ml of iodine. The range 2.5 -10 ml of 3% w/v Lugol's therefore corresponds to between 0.0632 and 0.2529 g of iodine, i.e. 0.4982-1.9928 mmol. The initial concentration of molecular iodine in Lugol's is unknown, so the quantity f ¼ b 0 /m 0 will be determined alongside the reaction rate via parameter fitting.
Vitamin C has a molar mass of 176.12 g mol 21 and hence a 1000 mg tablet diluted in 30-120 ml water yields stock solutions with a concentration of 0.18926-0.04731 mol l 21 . Hence 5 ml contains 0.9463-0.2366

Discussion
Clock reactions have long provided an instructive pedagogical example in chemistry education, in addition to their industrial and biochemical importance. This study focused on the vitamin C clock reaction, the dynamics of which are governed by substrate-depletion without autocatalysis. A fast vitamin C-dependent reaction converts iodine to iodide, depleting vitamin C in the process. A slow vitamin C-independent reaction converts iodide back to iodine. During a long induction period, the fast reaction dominates and very little molecular iodine is present in comparison to iodide. Once the vitamin C supply is exhausted, the slow reaction takes over and the molecular iodine level rises, which can be visualized by a colour change in the presence of starch. This system can be described effectively via techniques of matched asymptotic expansions, previously used successfully for the indirectly catalytic iodate-arsenous acid reaction [8] in addition to other systems. This analysis enables the construction of an approximate solution, and moreover a compact expression for the switchover time which depends on the initial concentrations of vitamin C and iodine, and the rate of the slow reaction, i.e. Results from 'kitchen sink chemistry' experiments were found to follow this model very closely, with the parameters k 0 and f being fitted simultaneously to data series varying in each of c 0 and m 0 . The fitted parameter representing initial molecular iodine proportion f was approximately zero. The analysis presented here should be of interest in both applications involving substrate-depletion dynamics, and also for pedagogical purposes in mathematical chemistry. The core ideas employed: law of mass action, dimensional analysis, quasi-steady kinetics, phase planes, matched asymptotics, numerical solutions and parameter fitting, are central to mathematical biology and chemistry, and within the reach of advanced undergraduates and masters students in mathematical modelling. Moreover, the experiment can be carried out using relatively safe chemicals and minimal equipment. A potential interesting avenue for further investigation, also accessible to student modellers, would be to vary the temperature of the reactants and assess whether the change to switchover time may be predicted via an Arrhenius equation for k 0 . A further topic to explore would be to use techniques from analytical chemistry such as UV-visible absorbance [19]   royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181367