Mathematical modelling of the vitamin C clock reaction

‘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 sulphite–iodate reaction in the 1880s – for review see Horváth & Nagypál [1], and the work of G. Forbes [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 multistep reactions via systems of nonlinear differential equations was given by Chien [6]; this work exhibited a number of analytical techniques including solving a Riccati-type 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 references [1] and [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\rightarrow 3B$, where $B$ is the clock chemical and $P$ is a precursor. This type of reaction is characterised 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\rightarrow D$, keeping the concentration of $B$ low. Simultaneously, clock chemical is produced upstream from a precursor, e.g. via the reaction $P\rightarrow 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, characterising the process as involving a combination of inhibition ($B+C\rightarrow 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\rightarrow 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₂), 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^-). Further mathematical development can be found in references [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 reference [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,

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^{-}$ and $H_{2}O_{2}$ concentrations due to a rate-limiting step involving nucleophilic attack of $I^{-}$ 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^{-}$ 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,

or in simplified form,

In contrast to the iodate-arsenous acid system, the pair of reactions (2), (4) are not overall autocatalytic in either $A$ or $B$; the sum $a+2b$ is conserved. The reaction has similar form to the reaction of iodide and peroxydisulphate in the presence of thiosulphate, which has been characterised by Horváth and Nagypál [1] as a pure substrate-depletive clock reaction (and therefore fitting Lente et al.’s strict definition of a clock reaction [13]). When sufficient inhibitor $C$ is present, the reaction (4) dominates, and iodine is held at relatively low concentration in comparison to iodide ions. However, the reaction (4) depletes vitamin C, and hence eventually this reaction can no longer continue. At this point, the slower reaction (2) dominates, resulting in concentration of the clock chemical iodine rising, and hence the solution turning from clear/white to blue. The time at which this reversal in the dynamics occurs is repeatable. A similar reaction involving persulphate in place of hydrogen peroxide has recently been analysed by Burgess and Davidson [19].

From the law of mass action, the reactions (2), (4) and the assumption that the system is well-mixed and isothermal, lead to the ordinary differential equation model,

with initial conditions based on the initial concentrations of $A$, $B$ and $C$ respectively.

Inspecting equations (5), (6) 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 2-variable system,

Nondimensionalising the system with the scalings,

leads to the dimensionless initial-value problem,

with initial conditions,

where the dimensionless parameters $\rho=m_{0}/c_{0}$ and $\epsilon=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. $\epsilon\ll 1$. The ratio $\rho$ will be assumed to be order $1$; it will also turn out to be that $\rho\phi:=b_{0}/c_{0}<1$ in order for the vitamin C supply to survive the initial transient.

Significant insight into the induction period of the system (13) (14) 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 $d\beta/d\tau\approx 0$ in equation (13), then it follows that $\beta\gamma\approx\epsilon\rho(1-2\beta)^{2}$. Substituting this expression into equation (14) for inhibitor depletion leads to,

Given that clock chemical concentration is small during the induction period, $\beta\approx 0$ and so equation (16) can be simplified as $d\gamma/d\tau\approx-\epsilon\rho^{2}$, therefore $\gamma\approx c_{1}-\epsilon\rho^{2}\tau$ where $c_{1}$ is a constant which we will determine in section 4.2. This expression shows that the length of the induction period scales with $(\epsilon\rho^{2})^{-1}$ in dimensionless variables.

The system (13), (14) has the unique equilibrium $(\beta,\gamma)=(1/2,0)$ which represents the long term fate of the system (zero inhibitor, all reactants converted to clock chemical). The eigenvalues of the linearised system at this point are $0$, with eigenvector $(1,0)$ and $-\rho/2$, with eigenvector $(1,\rho)$. The zero eigenvalue indicates the slow manifold $\{(s,0):s\in\mathbb{R}\}$ which we will find the system approaches after the induction period is complete. The negative eigenvalue indicates a stable manifold $\{(s,\rho s):s\in\mathbb{R}\}$ which lies outside the physically-reasonable region of interest $0\leqslant\beta\leqslant 1/2$, $0\leqslant\gamma$.

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

In this section we present the results of some ‘kitchen sink’ experiments, conducted with readily-available chemicals. The protocol is essentially as described by [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.

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 40^∘C. 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 for iodine and 166 g/mol for potassium iodide, Lugol’s 3% w/v solution contains $1+2\times 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 $\phi=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 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. Hence 5 ml contains 0.9463–0.2366 mmol of vitamin C. Concentrations are calculated based on 120 ml water, 5 ml stock, 15 ml peroxide and 2.5–10 ml Lugol’s. The concentration of hydrogen peroxide in all experiments is $0.098$ mol/l.

Results for two series, (a) varying vitamin C with iodine held fixed, and (b) varying iodine with vitamin C held fixed, are shown in table 1. The outcome of unconstrained least squares fitting Equation (29) for the parameters $k_{0}$ and $\phi$ to both experimental series simultaneously is shown in figure 5; we find $k_{0}\approx 0.57$ M^-1s^-1 and $\phi=-7\times 10^{-5}\approx 0$.

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 visualised 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 $\phi$ being fitted simultaneously to data series varying in each of $c_{0}$ and $m_{0}$. The fitted parameter representing initial molecular iodine proportion $\phi$ 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] to measure the trajectory of the solution quantitatively and compare to the mathematical solution. Clock reactions provide enduring and accessible examples of mathematical modelling and experiment.

Acknowledgements: DJS acknowledges Engineering and Physical Sciences Research Council award EP/N021096/1.

Data accessibility: Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.68q4hf7