Mathematical modelling of the vitamin C clock reaction: a study of two kinetic regimes
A. A. Alsaleh, D. J. Smith, S. Jabbari

TL;DR
This paper models a vitamin C clock reaction to better understand how hydrogen peroxide affects its timing and behavior.
Contribution
A more detailed and validated mechanistic model of the vitamin C clock reaction is developed and analyzed.
Findings
The model explains the effect of hydrogen peroxide concentration on the clock reaction's switchover time.
The model resolves discrepancies in the literature about the slow reaction's kinetic order.
The model is validated across a wide range of hydrogen peroxide concentrations.
Abstract
Chemically reacting systems exhibiting a repeatable delay period before a visible and sudden change are referred to as clock reactions; they have a long history in education and provide an idealization of various biochemical and industrial processes. We focus on a purely substrate-depletive clock reaction utilizing vitamin C, hydrogen peroxide, iodine and starch. Building on a recent study of a simplified two-reaction model under high hydrogen peroxide concentrations, we develop a more detailed model which breaks the slow reaction into two steps, one of which is rate-limiting unless hydrogen peroxide levels are very high. Through asymptotic analysis, this model enables the effect of hydrogen peroxide concentration to be elucidated in a principled way, resolving an apparent discrepancy with earlier literature regarding the order of the slow reaction kinetics. The model is analysed in…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5- —Kingdom of Saudi Arabia, Ministry of Education
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsChemistry and Chemical Engineering · Ionic liquids properties and applications · Molecular Junctions and Nanostructures
Introduction
Clock reactions are characterized by a defined and predictable induction period followed by a sudden and typically visible change in reactant concentrations. The study of these reactions dates back at least to the work of Landolt on the sulfite/iodate reaction in the 1880s [1]. These systems have long been used in chemistry education [2], have industrial applications [3,4] and alongside experiment have been studied through mathematical and computational models [5–9]. Clock reactions are of current interest through recent applications as diverse as a chemical clock car activity used in an educational setting [10], the evaluation of three-dimensional printed mixing devices [11] and the determination of microconcentrations of the potentially toxic dye indigo carmine [12]. For further details, see the reviews [13,14].
In this article, we will focus on a specific version of the vitamin C clock reaction, described by Wright [15,16] with the aim of providing a safe version of the experiment using only the ‘household chemicals’ of iodine, hydrogen peroxide, starch and ascorbic acid (vitamin C). The system involves two opposing reaction processes: a ‘fast reaction’ which consumes vitamin C ( ) and converts molecular iodine ( ) to iodide ions ( ),
and a ‘slow reaction’ which consumes hydrogen peroxide ( ) and converts iodide ions to molecular iodine,
As will be discussed below, the overall reaction (1.2) encapsulates the effect of a rate-determining reaction between iodide ions and hydrogen peroxide to produce hypoiodous acid, and the second faster step involving the combination of hypoiodous acid with a second iodide ion to produce molecular iodine [17].
Following initial mixing of iodine, hydrogen peroxide, starch and vitamin C, the fast reaction initially dominates, ensuring that the iodide concentration is much higher than the iodine concentration; in the presence of starch, the solution appears white. The solution remains in this state until the vitamin C concentration is nearly depleted, at which point the slow reaction then dominates, converting iodide ions to molecular iodine. In the presence of starch, the colour of the solution darkens so that it appears blue. The time interval at which this colour change occurs will be referred to as the ‘switchover time’ ; this interval is repeatable, meaning that the system fits Lente’s strict definition of a clock reaction [18]. Within the terminology of Horváth & Nagypál [13], the system is purely substrate depletive, by contrast with systems which involve autocatalysis or other clock mechanisms. The subject of this article is predicting on the basis of a mathematical model how depends on the initial concentrations of the reactants. For other clock reactions involving vitamin C, see references [11,19].
The vitamin C clock reaction has been modelled mathematically and tested experimentally [9]; this work focused on a regime in which hydrogen peroxide levels are greatly in excess compared with other reactants, which enabled a simplified model involving only two variables, iodine and vitamin C. By contrast, a recent pedagogical application [10] worked in a regime in which hydrogen peroxide levels were comparable to the other reactants. In addition to requiring the inclusion of an additional variable to take account of the consumption of hydrogen peroxide, this regime makes a qualitatitive change to the apparent kinetics of the slow reaction, from being effectively quadratic in the high hydrogen peroxide case to effectively linear with moderate hydrogen peroxide. The question of whether the reaction should be modelled as linear or quadratic arose in the peer review of reference [9], which proved a significant point of contention due to previous experimental work showing that this reaction is linear [17,20,21]. As discussed in [9] and the associated open peer review, quadratic reaction kinetics were found to lead to a formula for the switchover time which much better fitted the experimental data of [9] than a linear version. It was hypothesized that the reason for this discrepancy was the major difference in hydrogen peroxide concentration used by Kerr et al. [9] relative to these earlier papers. The hypothesis was that the hypoiodous acid generating step is rate-limiting only when hydrogen peroxide levels are not in great excess.
In this article, we will develop a mathematical model which accounts for the details of the production and conversion of hypoiodous acid within the slow reaction, thereby providing a unifying framework for both moderate and high hydrogen peroxide concentrations. By a similar approach to Billingham & Needham [7] and Kerr et al. [9], we will analyse the system through matched asymptotic expansions, exploiting the disparity of slow and fast reaction rates, to develop approximate analytical expressions for the timecourse of the reactant concentrations from initial mixing to final equilibrium. This analysis will lead to approximate formulae for the switchover time in moderate and high hydrogen peroxide regimes in terms of the initial concentrations, one reaction rate parameter, and a parameter relating to the initial iodide/iodine ratio (thus also identifying the dominant reactions and concentrations in the system that drive the switchover time). The mathematical models will be tested through their ability to fit experiment data with the same shared parameter values.
Model formulation
As in [9], the fast reaction (1.1) will be expressed in terms of model variables as,
where is concentration of molecular iodine, is concentration of vitamin C, is concentration of iodide and is the associated reaction rate, as in [9].
The departure point for the present work is in separating the slow reaction (1.2) into two steps,
where the first step is the production of hypoiodous acid , and the second step completes the overall conversion of two iodide ions to one molecule of iodine. The first step is taken to be rate-limiting under conditions of moderate hydrogen peroxide concentration, which can be expressed mathematically as .
Introducing the variables for hydrogen peroxide and for hypoiodous acid, these reactions are simplified as,
(note that is not considered in the model as it has no influence on the other variables).
Finally, there exists a reverse pathway by which hypoiodous acid is reduced by hydrogen peroxide to yield iodide [22],
or in model variables,
Applying the law of mass action to equations (2.1), (2.3a), (2.3b) and (2.5) yields the system of ordinary differential equations,
with initial conditions expressed as
In summary, the model contains five dependent variables and has four free parameters for the reaction rates. The initial values and are straightforward to control experimentally, as is the total atomic iodine ; it will be assumed below that the initial concentration of hypoiodous acid . We will therefore consider to be a fifth free parameter.
Note that conservation of atomic iodine, is constant and equal to its initial value . The model can therefore be reduced to
We will make the following assumptions regarding the relative rates of the reactions in the system. The fast reaction is assumed to take place one order of magnitude faster than the fastest part of the slow reaction. The rate-limiting step (2.3a) is assumed to be one further order of magnitude slower, as is the reverse reaction (2.5); schematically,
More precisely, nondimensionalizing with scalings
where denotes a dimensionless quantity, yields the dimensionless system,
with initial conditions
The dimensionless parameters are then defined as,
where , , are order 1, and quantifies the disparity of the reaction rates. The magnitude of will define moderate (order 1) or high (order ) hydrogen peroxide regimes. Note that the definition of is slightly different from that used in [9], which expressed the rates in terms of the rate constant of a simplified form of the slow reaction. In what follows, asterisks will be dropped for brevity.
As in earlier studies [7–9], we will exploit the small parameter to develop an approximate solution to the problem via matched asymptotic expansions. The regimes of moderate and high hydrogen peroxide will be denoted M-HP and H-HP, and specifically refer to whether the parameter is order 1 or order , which determines whether the hypoiodous acid-producing step has magnitude or in the dimensionless system. We will now proceed to analyse each regime region-by-region. A comparison of numerical solutions to equation (2.10) (computed with Matlab ode23s) and the asymptotic solutions which will be constructed is given in figure 1. The switchover time can be seen at the point where vitamin C (red line) drops, inducing a sharp increase in iodine (blue line).
Comparison of numerical solutions (solid lines) and asymptotic solutions (dotted lines), as found in §3 and §4, respectively, to the dimensionless problem (2.10) plotted against logarithmic time. Parameter values (chosen arbitrarily): ϵ=10−2 , β=0.6 , γ=0.7 , σ=0.8 and ϕ=0.2 . (a) Moderate hydrogen peroxide case M-HP with ρ=2 ; (b) high hydrogen peroxide case H-HP with ρ^=0.9 so that ρ=ρ^ϵ−1=90 .
Model M-HP: moderate hydrogen peroxide
ρ=0(1)
Under the assumption that is order 1 (i.e. hydrogen peroxide concentration is initially comparable with vitamin C concentration), we note that hypoiodous acid is initially zero and its rate of production is order . Therefore, it is appropriate to rescale where is order 1. Then the rescaled system is,
with initial conditions (2.11). The Model M-HP (3.1) will be studied through four asymptotic regions in detail in the following sections.
Region I: initial adjustment
3.1.
During region I, the reactants are mixed before the induction phase begins. It can be described as where the independent variable is order 1 ( ) and the dependent variables are order 1. Seeking a solution in the form of asymptotic expansions:
the leading order terms for the asymptotic expansion of the Model M-HP (3.1) are:
with the initial conditions:
Equations (3.2a) and (3.2b), alongside their initial conditions, yield
Dividing (3.2c) by (3.2d) gives
which has solution
Substituting into equation (3.2d) then gives the differential equation in terms of only,
which has solution,
Substituting (3.6) into (3.5) leads to
Throughout, it will be assumed that in dimensional variables, , which corresponds to the nondimensional condition . This condition corresponds to the initial vitamin C concentration being sufficient to survive the initial adjustment. The initial adjustment therefore corresponds to consumption of a proportion of the initial vitamin C concentration via the fast reaction, thereby converting the bulk of the initial iodine concentration to iodide in order 1 time; the slow reaction is subleading in this region. In particular, note that the order 1 part of is exponentially decaying.
Region Ia: quasi-equilibrium established for
Q
3.2.
Following the exponential decay of vitamin C (to ) on the previous region, the fast reaction no longer dominates and the slow reaction becomes of comparable importance. By inspecting equation (3.1d), it is apparent that balance between fast and slow reactions implies that .
A distinctive feature of the M-HP case as opposed to the H-HP case is the existence of a region intermediate between the initial adjustment and the induction period, during which the hypoiodous acid concentration reaches quasi-equilibrium following the rise in iodide concentration. For consistency with [9], which used the designation ‘region II’ specifically for the induction period, we will designate the intermediate region between the initial adjustment and the induction period by ‘Ia’. Physically, region Ia corresponds to the two forward steps in the slow reaction evolving towards a balance in their production and removal of hypoiodous acid. Inspecting equation (3.1b) shows that varies on a timescale . Therefore, the scalings for region Ia are, , and . Denoting and , then the system takes the form
The leading order terms of the asymptotic expansions therefore satisfy,
with the matching conditions to the region I solution:
Solving equations (3.9a) and (3.9c) then yields the constant solutions,
Equations (3.9b) and (3.9d) can then be solved to yield,
and
The limiting behaviours of the solution in region Ia are therefore,
as .
Region II: induction period
3.3.
Following the establishment of quasi-equilibrium within the slow reaction, the system enters the induction period, which is the distinguishing feature of a clock reaction. During the induction period, the fast and slow reactions are in quasi-equilibrium, with iodine levels held low due to the disparity in reaction rates. The vitamin C is gradually consumed until it is exhausted and the switchover process then occurs. From inspecting equation (3.8c), it is apparent that this process occurs on a timescale of order slower than , motivating the introduction of a second rescaled time variable . The rescaled system is then,
The leading order terms of the asymptotic expansions in region II, therefore, satisfy,
Solving equation (3.15) and matching to the region Ia solutions then yields the leading order approximations,
The switchover time can then be approximated as the time at which the region II solution predicts that the vitamin C concentration is zero at leading order. Inspecting equation (3.16c) yields in dimensionless variables,
or in dimensional variables,
This equation emphasizes that there must be sufficient combined hydrogen peroxide and molecular iodine provided initially to drive the vitamin C concentration to zero before the hydrogen peroxide is exhausted, i.e. .
A similar equation to (3.18) was given in [10] through a quasi-steady argument and making the phenomenological assumption that the slow reaction is linear in iodide concentration (although not including the term in the denominator). However, by accounting for the sub-steps of the slow reaction, the dynamics of this regime emerge naturally from the law of mass action. Equation (3.18) is one of the two formulae that will be tested experimentally in this paper.
While the above analysis provides a formula for the switchover time, as in [9], it is of interest to construct approximate solutions for the remaining dynamics of the system.
Region III: corner
3.4.
As the vitamin C concentration is exhausted, the iodine concentration begins to change significantly, and the system is no longer in equilibrium. Following [9], we will refer to this as the ‘corner’ region; the analysis is similar, with some minor adaptations for the four variable model and different definition of the small parameter .
Shifting the time coordinate so that the origin corresponds to the switchover time, the most structured balance in equation (3.1) is given by the ‘crossover’ scalings,
The rescaled system is then,
Taking the leading order terms of the asymptotic expansions yields,
Matching (3.21a) to region II at gives . Substituting into equation (3.21b) and solving for , we arrive at
where is a constant. Matching to region II as establishes that so that remains finite; indeed the resulting constant solution then matches with the region II solution at .
Subtracting a multiple of equation (3.21d) from (3.21c) and substituting the solution for then gives the equation,
from which it follows that
where is a constant. Noting that an order 1 change to corresponds to an order 1 shift in the origin for , and hence a subleading change to , it is reasonable to choose for convenience in the present analysis. A more precise choice of would require the next order solution for in region II to be calculated for the purpose of matching. Using equation (3.24) to substitute for in equation (3.21c) then yields a Riccati equation for ,
as found for example in the study of a cubic autocatalysis clock reaction [7]. It is convenient to introduce the notation which renders the equation into an identical form to that in [9], from which we can immediately write down the solution,
where is another constant to be determined. The constant can be deduced by examining the limiting behaviour of the error function as ; for the solution for to approach a non-zero limit requires that . The solution for follows immediately from equation (3.24). In summary, the solutions for region III are,
Region IV: equilibration
3.5.
Following the exhaustion of vitamin C and the switchover, the system approaches its long-term equilibrium state. From equation (2.10), it is clear that equilibrium with and also requires that and — in other words, all iodide and hypoiodous acid is ultimately converted to molecular iodine.
Using the same timescale as the induction period, but with now order 1, the rescaled system is,
Equation (3.28c) yields that at leading order, from which we deduce the expected result that for all . Indeed it can be shown by mathematical induction that at all powers of , by a similar argument to that given by Kerr et al. [9]. The leading order terms of the rest of the system are,
Rearranging this system yields,
where is a constant. Matching to region II immediately yields , i.e. the level of hydrogen peroxide remaining after the induction period. We will retain the notation for the rest of the calculation for brevity. Substituting for in equation (3.29a) then yields an equation involving only ,
which integrates to,
where is another constant to be determined by matching. Noting that and moreover at the beginning of region IV, , it follows that the term inside the modulus signs is initially positive; a change of sign of the numerator cannot then occur at finite time, so we deduce that,
The constant must then be,
Rearranging equation (3.33) then gives the long-timescale solution for hydrogen peroxide,
and substitution into equation (3.30) shows moreover that,
Two cases can be distinguished in equation (3.35). If , i.e. there is sufficient hydrogen peroxide remaining following the switchover, then in the limit as ,
in other words, all of the iodide and hypoiodous acid are converted to molecular iodine.
If however , i.e. there is insufficient hydrogen peroxide remaining following the switchover, then as ,
In dimensional variables, the condition is equivalent to . To convert all iodide to iodine by the end of the process, there must be sufficient hydrogen peroxide to account for all supplied atomic iodine (noting one hydrogen peroxide molecule is required for every two iodide ions), plus exhausting all supplied vitamin C, minus the initial adjustment.
In summary, the region IV leading order solutions describing the long-term evolution of the system are,
where is given in equation (3.17) and .
A comparison between numerical and leading order asymptotic solutions in each region in dimensionless variables is given in figure 1a, for a specific set of parameters, with small parameter chosen as . In the absence of detailed measurement of the reaction rates in the system, this choice is arbitrary. Additional detail on the dynamics for in region Ia is given in appendix A, figure 6. As shown in figure 1a, there is a small but noticeable difference between the asymptotic and numerical solutions as approaches zero in the induction region which results in a small discrepancy in switchover time (figure 2a). By plotting as a function of , the relative difference in the asymptotic switchover time value equation (3.17) and the point at which the numerical solution falls below a threshold value of (figure 2b), this discrepancy converges to zero. As would be expected from taking only the leading order terms of the asymptotic expansion, the convergence is approximately linear in .
Comparison of asymptotic-derived formula and numerical solution in determination of the dimensionless switchover time for regime M-HP. Parameter values (chosen arbitrarily): ϵ=10−2 , β=0.6 , γ=0.7 , σ=0.8 , ϕ=0.2 and ρ=2 . (a) Numerical solutions for C(t) and I(t) shown alongside the value of the switchover time (dashed line) calculated by the formula (3.17). (b) Relative error between the asymptotic expression (3.17) and the point at which the numerical solution falls below a threshold value of ϵ , as a function of the small parameter ϵ .
Model H-HP: high hydrogen peroxide
ρ=O(ϵ−1)
The contrasting regime with high hydrogen peroxide can be characterized mathematically by the ratio , or more precisely, being order . Therefore, we will write , where is order 1. The dimensionless model of equation (2.10) then becomes,
with the same initial conditions (2.11) as previously.
Region I: initial adjustment
4.1.
The initial process of vitamin C being rapidly consumed to convert the bulk of the iodide to molecular iodine occurs during a time interval of order 1 similarly to the M-HP regime. Again, it is necessary to rescale the hypoiodous acid concentration, so that , with being order 1. Then the rescaled system is,
The leading order system becomes,
The solutions for , and are identical to those in §3.1. The problem for then reduces to the separable equation,
with solution,
The limiting behaviour as is then,
The linear growth of implies that becomes order 1 on a timescale of order .
The more rapid growth of , arising because the hypoiodous acid-producing step is no longer rate-limiting in case of abundant hydrogen peroxide, means that it is no longer necessary to consider a separate region for quasi-equilibrium to be established, so there is no counterpart to region Ia studied for the M-HP regime.
Region II: induction period
4.2.
Moving on to region II—the induction period—it is apparent that quasi-equilibrium for forward and backward reactions occurs in equation (4.1d) when is order , i.e. an order of magnitude larger than in the moderate hydrogen peroxide case. Considering equation (4.1c) with the rescaling (with being order 1), it is clear that the appropriate time rescaling is then , with being order 1. The resulting system for region II is then,
The leading order system takes the form,
After integrating equation (4.8a) and matching to region I, the hydrogen peroxide concentration remains constant at leading order, i.e. . Substituting into equation (4.8b) and solving yields,
where the constants and are given by,
and is to be determined via matching. Because , it can be deduced that is real-valued.
Matching to region I (as , ), gives
from which we deduce the solution,
Combining equations (4.8c), (4.8d) and substituting the solution for then reduces the problem for to
which integrates to yield,
The constants are given by,
and the constant of integration is found by matching to the region I asymptote , giving the value .
Substituting equations (4.12), (4.14) into (4.8d) and rearranging for then completes the leading order solution in region II. This expression is quite lengthy to write out in full so the solution is expressed below in terms of and . The leading order solutions in region II are therefore,
The switchover time can again be approximated by solving . While one could in principle solve this equation numerically for , we will instead carry out a series of approximations to yield a closed-form expression. First, we assume that is sufficiently large to make the approximations and , which leads to the equation,
Therefore,
Equation (4.18) follows from a controlled approximation applied to the leading order solution, and could in principle be used for fitting experimental data. However, we will instead attempt to simplify this expression further to reduce the number of free parameters, and also to make contact with the result derived by Kerr et al. [9] under the assumptions of quadratic kinetics for the slow reaction.
The numerator of equation (4.18) can be written in terms of the dimensionless parameter groupings and as,
As mentioned above, , therefore the logarithmic and square root terms may be expanded in powers of , which gives the following leading order approximation for the numerator in equation (4.18),
We will now consider conditions under which the term in square brackets in equation (4.20) can be neglected. Considering all positive values for , the leading order term is bounded by its value when , in which case the term takes the value . Therefore, provided that is small in comparison with , the numerator of equation (4.18) can then be approximated by . In dimensional parameters, this assumption corresponds to , and so will be valid provided that the initial vitamin C concentration is chosen sufficiently high.
Under this assumption, the dimensionless expression then simplifies to,
or in dimensional variables,
Equation (4.22) can be considered a somewhat simplified expression for the switchover time, the simplifications being transparently (if not unconditionally) justified. Because there are four free parameters ( ) we will find it convenient to make one further assumption, which is that the denominator is dominated by the term. This assumption is uncontrolled and cannot be justified a priori, but has the considerable advantage of reducing the expression to the simple form,
Equation (4.23) has a similar structure to that of Kerr et al. [9] and like equation (3.18) only possesses two fitting parameters, and . The reasonableness of this step will be demonstrated through simultaneous fitting to experimental data in §5.
Because equation (4.22) includes an inverse dependence on , the question may be posed (see open peer review accompanying this article) whether this dependence predicts that switchover time can be made arbitrarily small in practice simply by increasing . However, sufficiently large values of would correspond to no longer being merely , which would be outside of the range of validity of the present analysis by changing the relative orders of magnitudes of the reaction steps. The implications of taking at the next order of magnitude are discussed briefly in appendix B; in short, the analysis predicts that the switchover time approaches a non-zero value which has no quantitative dependence on .
Region III: corner
4.3.
To connect the induction period to the long-term state it is again necessary to consider a corner region centred at . The most structured balance of the high peroxide system occurs with timescale and dependent variables scaled as and , with and remaining order 1. The scaled system is then,
Expanding in powers of , the leading order system takes the form,
The variables and are therefore constant in region III, their values being determined by matching to region II. As in region II, . From inspection of equation (4.12) and the properties of the tanh function, it is clear that the region II solution for rapidly approaches the steady state value,
For brevity, we will define the remaining equations are,
As in §3.4, the variable can be eliminated to yield a Ricatti equation,
where is a constant of integration, which is again indeterminate at leading order, so it is taken to be zero. This equation takes the same form as that given in §3.4 and therefore the solution has a similar form, and can then be deduced from the relation .
The leading order solutions in region III are therefore,
with as defined in equation (4.26) et seq.
Region IV: equilibration
4.4.
As for the M-HP case, following the switchover, is again at most order and , are order 1. By contrast with the M-HP case, is also order 1, and the time variable is order . Setting and , the system describing the long-term dynamics after the switchover is then,
The leading order system is,
Again, it is clear that and that ; indeed for all . The dynamics of and are therefore described by the pair of equations,
We have not been able to find a closed-form solution of the leading order equations (4.32a)–(4.32b) for order 1 values of the dimensionless parameters. The system has a single equilibrium point , which physically corresponds to all iodide and hypoiodous acid having been converted to iodine. The eigenvalues of this equilibrium point are and . The zero eigenvalue corresponds to the slow manifold which is tangent to the line . The negative eigenvalue corresponds to a stable manifold tangent to which is outside of the physical domain of interest because it would entail either or being negative.
The ‘terminal dynamics’ close to the equilibrium point can be inferred approximately by considering a quasi-steady balance of the terms linear in and in equation (4.32a), yielding,
which can be rearranged to provide the following approximate relationship between and ,
from which we deduce that,
Substituting equations (4.33) and (4.35) into (4.32b) then provides an approximation to the long-term dynamics of :
This equation has solution,
where is a constant of integration. The corresponding approximation to then follows from (4.34),
The constant can be determined by matching the solution for to the region II solution, giving the value . The approximate solutions in region IV, expected to be increasingly accurate close to the equilibrium point, are therefore,
Figure 1b compares the high hydrogen peroxide asymptotic solutions in each region with a numerical solution computed with an arbitrary set of parameters (details in caption). The solutions are quite close in their region of validity, particularly for the trajectories of , and , although as time progresses, a significant drift develops between the numerically computed value of and the constant approximation produced by the asymptotic analysis; the discrepancy appears to be on the order of . This discrepancy could be reduced by seeking the next order term in the asymptotic expansion, although this quantity is less critical to determine than the clock chemicals and which directly determine the switchover time.
Figure 3a compares the asymptotic-derived switchover time formula (4.18) to that of the numerical solution with the relative error between these for varying given in figure 3b. By contrast with the moderate hydrogen peroxide case, the convergence of the error as is sublinear, the gradient on a log–log plot approaching approximately . These results provide evidence that the asymptotic expansion has sublinear accuracy in this region and so a correction with a fractional power in may be needed. We will not pursue this analysis in the present paper as it appears likely to involve significant additional complications for limited additional insight.
Comparison of asymptotic-derived formula and numerical solution in determination of the dimensionless switchover time for regime H-HP. Parameter values (chosen arbitrarily): ϵ=10−2 , β=0.6 , γ=0.7 , σ=0.8 , ϕ=0.2 and ρ^=0.9 . (a) Numerical solutions for C(t) and I(t) shown alongside the value of the switchover time (dashed line) calculated by the approximate formula (4.18). (b) Relative error between the asymptotic expression for the switchover time (4.18) and the point at which the numerical solution falls below a threshold value of ϵ , as a function of the small parameter ϵ .
Overall, the comparison for each of the M-HP and H-HP problems in figure 1 provide strong evidence for the validity of the asymptotic approach. In particular, the approximations for the clock chemical iodine are both very satisfactory, as are the results for the inhibitor vitamin C, and the intermediate species hypoiodous acid. For large time (region IV for M-HP and late in region II for H-HP), a discrepancy emerges between the numerical solution and leading order asymptotic solution for hydrogen peroxide, likely due to the need for more terms or possibly a slower-decaying asymptotic series approximation. Because the error in hydrogen peroxide does not appear to affect the approximation to the clock chemical concentration or switchover time, we will not investigate this issue further in the present manuscript.
Experimental testing
The M-HP and H-HP switchover time formulae (3.18) and (4.23) were tested through tabletop experiments in a similar manner to the work of Kerr et al. [9], using a combination of vitamin C powder, Lugol’s iodine, hydrogen peroxide solution, and powdered laundry starch. A key refinement to the experimental technique was to use a webcam sensor running under Matlab R2020b to detect the point at which the colour change occurred. All experiments were carried out diluted in 60 ml water at 40°, with 5 ml of vitamin C stock solution (made up as 1000 mg in 30−90 ml water), 3−10 ml of 3% Lugol’s iodine, 5 g laundry starch and 1−20 ml 3% hydrogen peroxide.
Imaging
5.1.
The switchover time was measured by imaging the mixture from above under natural lighting conditions using an RGB USB camera running under Matlab with Image Processing Toolbox (Mathworks, Natick) 2020b, recording images at 15 frames per second. The region of interest was set at pixels in the centre of the beaker. The red channel was sufficient for a clear measurement. Briefly, the signal was processed by summing the intensity across all pixels, taking the difference of successive frames, then taking forward and backward moving averages with a 10-frame window to minimize the effect of fluctuations, then identifying the point at which the difference between forward and backward moving averages took the largest negative value. This ‘corner’ value corresponds approximately to the point of largest second derivative of the signal, and is therefore relatively insensitive to lighting conditions. Image processing code can be found in the accompanying supplement [23].
Fitting series
5.2.
Four experimental series, referred to as and (moderate hydrogen peroxide regime), and and (high hydrogen peroxide regime) were carried out for fitting of the constants and . These experiments involved holding varying respectively initial iodine concentration or vitamin C concentration while holding the initial masses of other substances fixed. Due to the slightly varying volume of the overall solution, in series and , the concentrations of vitamin C and hydrogen peroxide also varied slightly, which was taken into account when fitting. Full details of the concentrations aare given in the accompanying supplement [23].
The moderate hydrogen peroxide regime involved using ml hydrogen peroxide, providing an initial concentration of mol l^−1^. In series , the initial mass of Lugol’s was varied from to ml, producing an initial concentration of to mol l^−1^, with mol l^−1^. In series , the initial mass of vitamin C was varied from to g, producing a concentration of to mol l^−1^, with mol l^−1^. In all of these experiments, the initial ratio of hydrogen peroxide to iodine was in the range . Measured values of the switchover time varied between approximately 60 and 700 s.
The high hydrogen peroxide regime used ml hydrogen peroxide, providing an initial concentration of mol l^−1^. In series , the initial mass of Lugol’s was varied from to ml, producing an initial concentration of to mol l^−1^, with mol l^−1^. In series , the initial mass of vitamin C was varied from to g, producing a concentration of to mol l^−1^, with mol l^−1^. In these experiments the ratio . Measured values of the switchover time varied from 23 to around 190 s.
The values of and were estimated through relative error least squares fitting in order to account for the wide range of measured switchover time, carried out in Matlab with the unconstrained local optimization function fminsearch, which implements the Nelder–Mead simplex search method [24]. Confidence limits were estimated through bootstrapping with samples, implemented in Matlab with the function bootci, which implements the bias corrected and accelerated percentile method [25].
Testing series
5.3.
The switchover time model and fitted parameters were then tested by comparing against two independent data series, referred to as and . Series varied the hydrogen peroxide concentration between and mol l^−1^, while holding and mol l^−1^, so that . Series varied the hydrogen peroxide concentration between and mol l^−1^, while holding and mol l^−1^, so that . The switchover time formula (3.18) was used for cases where and (4.23) when , the threshold value of being chosen arbitrarily.
Results
5.4.
Figure 4a–d shows the outcome of the fitting series experiments , , and . Qualitatively, the results show that in both low and high hydrogen peroxide regimes, switchover time declines as initial iodine concentration is increased, and increases nearly linearly with initial vitamin C concentration—fundamentally, as substrate is increased, inhibitor is used up more quickly, whereas as inhibitor is increased, it lasts longer. The quality of the fit appears excellent, with the maximum relative error rarely exceeding 10%. The estimated parameter values of (95%confidence interval ) and (95% CI ) are shown in figure 4e along with bootstrap sample distributions.
Outcome of simultaneous fitting to data series NM , CM , NH and CH . (a–d) experimental (blue dot) and fitted switchover time (red line) with sum square relative error best fit of k^2=0.066302 and ϕ^=0.15835 . (e) Bootstrapping results with 10 000 repeats (blue dots), best fit (red circle) and 95% confidence intervals for each parameter (dashed red lines).
The testing series and using the estimated parameters provide at least as close a match between experiment and the fitted model (figure 5) as for the fitted data, providing confidence in the out-of-sample applicability of the model and accuracy of the fitted parameters. Qualitatively, the data show switchover time decreasing as hydrogen peroxide concentration is increased, corresponding to the slow reaction being promoted through additional availability of a key reactant in its rate-limiting step.
Test of the parameter fits (figure 4) against independent experimental series H1 and H2 in which initial hydrogen peroxide mass is varied.
Discussion
This paper unifies and significantly extends previous work on asymptotic analysis and experimental data fitting for a tabletop chemical kinetics experiment. Through breaking down the reaction converting iodide to iodine into two steps involving the intermediate hypoiodous acid, the resulting model was able to explain the kinetics observed in regimes of moderate and high hydrogen peroxide concentration, encompassing both the models contained in Kerr et al. [9] and Parra Cordova & González Peña [10], and being transparently consistent with previous descriptions of the kinetics of the iodide–hydrogen peroxide reaction [17,20,21]. The key kinetic insight provided by the model is that the production of hypoiodous acid is only rate-limiting in the situation where hydrogen peroxide is comparable with other reactants.
The analysis utilized a single small parameter which was used to represent the disparity of rates between the slow rate of production of hypoiodous acid from iodide, the production of iodine from iodide and hypoiodous acid, and the conversion of iodine back to iodide through oxidation of vitamin C (fastest). A slow reverse reaction reducing hypoiodous acid back to iodide was also included. The small parameter was thereby able to distinguish the moderate from high hydrogen peroxide regime, the hydrogen peroxide-to-vitamin C ratio being order 1 in the moderate case, and order in the high case.
Asymptotic analysis was carried out in each of the two regimes. The moderate hydrogen peroxide case involved five distinct time regions: an timescale (region I) during which the initial free iodine is mostly converted to iodide through reduction of a fraction of the initial vitamin C. This region is followed by an process (region Ia) during which iodide, iodine and hypoiodous acid reach quasi-equilibrium. When comparing against numerical solutions, the leading order terms captured the kinetics of these species well, although there was some divergence from the hydrogen peroxide and vitamin C concentrations in this region. The induction period (region II) followed, which occurred on time. This process captured the consumption of hydrogen peroxide and vitamin C quite accurately, and provided a closed-form solution (3.18) for the switchover time, through solving for the leading order solution for the vitamin C concentration reaching zero. This formula can be seen as a refinement to that given by [10] by taking into account the effect of initially non-zero molecular iodine in oxidizing vitamin C and hence shortening the induction period. By similar methods to Kerr et al. [9], it was also possible to construct an approximate solution in the intermediate region (III) for order time around the switchover, and the long-term evolution to equilibrium (region IV) on a timescale of order . The leading order solutions on regions II–IV showed reasonably good agreement with numerical solutions, in particular, the relative error between the switchover time formula and the point at which the numerical solution fell below was found to converge approximately linearly as . An area for possible future development would be to predict better the consumption of vitamin C and hydrogen peroxide in region Ia, and to provide a higher-order approximation in region II, however as noted above, it seems unlikely that significant additional insight would be obtained.
The second asymptotic analysis of the paper concerned the high hydrogen peroxide case already considered by Kerr et al. [9]. Heuristically one can argue that when hydrogen peroxide levels are sufficiently high, the rate-limiting step is no longer such, and so the overall reaction can be modelled by the quadratic kinetics that would follow from the law of mass action, and hydrogen peroxide levels need not be modelled explicitly as they are not significantly reduced. The analysis of §4 put this informal argument on a more quantitative basis through elucidating how a hydrogen peroxide concentration of order provides the balance described above. Similar to Kerr et al. [9], we found four asymptotic regions corresponding to (I) initial reduction of vitamin C due to non-zero initial iodine, which is simultaneous with iodide, iodine and hypoiodous acid reaching quasi-equilibrium, (II) induction, (III) crossover and (IV) evolution to final equilibrium. The leading order solution for region IV could not be solved exactly, however a linearization about the equilibrium combined with matching to region II enabled a relatively accurate approximate solution to be found; however the leading order analysis was not even in principle able to capture the significant consumption of hydrogen peroxide seen in region IV.
One of the goals of this paper has been to show how previously reported formulae for the switchover time can be recovered from our unified model, when considered in the appropriate concentration regime. As previously, the region II solution enabled the deduction of a ‘full’ switchover time formula (4.18) for the full H-HP regime. This full formula was shown to be equivalent to a simplified version that could be justified rigorously under certain assumptions regarding parameters (4.22) and a highly simplified heuristic version for comparison with Kerr et al. [9] (4.23). Due to involving only two free parameters, it was the latter formula that we subsequently used for fitting experimental data. The full formula was found to converge to the numerical solution as sublinearly, suggesting that the next order correction may be sublinear.
The switchover time formulae were then tested through fitting the parameters (dimensionless initial proportion of molecular iodine) and (production of hypoiodous acid from iodide and hydrogen peroxide). Parameters were fitted through minimising least square relative error in order to provide similar weighting to results across more than an order of magnitude in time values. Simultaneous fitting was carried out by varying the initial mass of iodine or vitamin C, repeated in moderate and high hydrogen peroxide regimes—four series to estimate two parameters. Uncertainties were quantified through bootstrapping, with the 95% confidence interval for each parameter being on the order of %of the estimated value. The formula and estimates were then assessed against two separate experimental series in which initial mass of hydrogen peroxide was varied from moderate-to-high values. The model curves fell very close to the experimental results without any further fitting or tuning, giving high confidence in both the model and estimates.
The paper thereby presents a validated mechanistic understanding of a paradigm chemical kinetics system. We believe the work may form a useful basis to develop quantitive understanding of other substrate-depletive clock reactions and their applications, in addition to providing an accessible case study in the combination of time-dependent asymptotic analysis and mathematical model parameterization, an approach which is increasingly used to study biological and biochemical systems [26]. There are several areas in which the study may be developed and improved further.
In respect of the ‘mathematical’ aspects, in matching the induction region to the corner region, it was noted that more accuracy could have been obtained if a higher-order approximation to the induction region solution were available; this issue may also underlie the sublinear convergence of the numerical and leading order solutions for the switchover time observed in the high hydrogen peroxide model. It was also noticeable that the leading order constant solution for hydrogen peroxide loses accuracy at large time values; it would be of interest to explore whether these dynamics could easily be captured through a higher-order asymptotic solution.
Focusing on data fitting, while the strategy taken in this paper was to carry out data fitting on a highly simplified form of the high hydrogen peroxide switchover time formula, it should in principle be possible to fit the full form following from equation (4.18). The fitted model parameters may also provide post hoc justification of the simplifications carried out to form equation (4.22).
From a chemical kinetics modelling perspective, the hypoiodous acid-producing reaction has been found to comprise two parallel pathways, one of which depends on the availability of ions [17]. Therefore, it would be of interest to model these pathways independently in order to predict the effect of pH on the overall reaction rate. While the present study fixed the temperature of the reaction, another area of development would be to model temperature dependence of reaction rates and fit the model parameters by varying temperature systematically.
In respect of experimental technique, another area for development would be the use of ultraviolet light spectrophotometry [27] to measure the time evolution of iodine levels (via blue iodine–starch complex) rather than simply the switchover time, so that model curves could be fitted directly, providing a more complete test of the accuracy of the solutions in each asymptotic region. Explicit mathematical modelling of the formation of starch–iodine complex may in turn be required to interpret spectrophotometry data. As the discussion above illustrates, clock reactions continue to provide challenge and inspiration at the interface of mathematics and physical chemistry.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Landolt H . 1887 Ueber die Zeitdauer der Reaction zwischen Jodsäure und schwefliger Säure. Berichte Der Dtsch. Chem. Ges. 20 , 745–760. (10.1002/cber.188702001173) · doi ↗
- 2Conway WJ . 1940 A modified lecture experiment. J. Chem. Educ. 17 , 398. (10.1021/ed 017p 398) · doi ↗
- 3Billingham J , Coveney PV . 1993 Simple chemical clock reactions: application to cement hydration. J. Chem. Soc. Faraday Trans. 89 , 3021. (10.1039/ft 9938903021) · doi ↗
- 4Billingham J , Coveney PV . 1994 Kinetics of self-replicating micelles. J. Chem. Soc. Faraday Trans. 90 , 1953. (10.1039/ft 9949001953) · doi ↗
- 5Chien JY . 1948 Kinetic Analysis of Irreversible Consecutive Reactions. J. Am. Chem. Soc. 70 , 2256–2261. (10.1021/ja 01186 a 078) · doi ↗
- 6Anderson JM . 1976 Computer simulation in chemical kinetics. J. Chem. Educ. 53 , 561. (10.1021/ed 053p 561) · doi ↗
- 7Billingham J , Needham D . 1992 Mathematical Modelling of Chemical Clock Reactions. I. Induction, Inhibition and the Iodate-Arsenous-Acid Reaction. Philos. Trans. R. Soc. Lond. Ser. 340 , 569–591.
- 8Billingham J , Needham DJ . 1993 Mathematical modelling of chemical clock reactions. J. Eng. Math. 27 , 113–145. (10.1007/bf 00127478) · doi ↗
