EASY TO USE PROGRAM “SIMKINE3” FOR SIMULATING KINETIC PROFILES OF MULTI-STEP CHEMICAL SYSTEMS AND OPTIMISATION OF PREDICTABLE RATE COEFFICIENTS THEREIN

‘Simkine3’, a Delphi based software is developed to simulate the kinetic schemes of complex reaction mechanisms involving multiple sequential and competitive elementary steps for homogeneous and heterogeneous chemical reactions. Simkine3 is designed to translate the user specified mechanism into chemical first-order differential equations (ODEs) and optimise the estimated rate constants in such a way that simulated curves match the experimental kinetic profiles. TLSoda which uses backward differentiation method is utilised to solve resulting ODEs and Downhill Simplex method is used to optimise the estimated rate constants in a robotic way. An online help file is developed using HelpScrible Demo to guide the users of Simkine3. The versatility of the software is demonstrated by simulating the complex reaction between methylene violet and acidic bromate, a reaction which exhibits complex nonlinear kinetics. The new software is validated after testing it on a 19-step intricate mechanism involving 15 different species. The kinetic profiles of multiple simulated curves, illustrating the effect of initial concentrations of bromate, and bromide were matched with the corresponding experimental curves.


INTRODUCTION
Chemical kinetics is concerned with the dynamics of temporal concentrations of reactants and products [1].Practical experiments involve determination of an average rate of reaction of large number of molecules for the rate-determining step [2].The rates of chemical reactions are function of the extent to which the reaction has proceeded and the rate is a function of the reactant and/or product concentrations; it is found in all cases to be directly proportional to the amounts of one or more reactants and/or products raised to some power, which is the order with respect to the species and could be established experimentally [3].The rate laws are differential equations because the rate of reaction is the rate of changes of the extent of the reaction with time and as a consequence it is necessary to use the medium of calculus to write and solve these equations [4].
Detailed chemical kinetic models incorporate elementary reactions coupled with transport of matter and heat.The models are suitable for many purposes and can be used to describe reactive flows such as those encountered in the industrial chemical processes or to better understand combustion or environmental chemistry [5].The elucidation of reaction mechanism involves a mix of experimentation, literature search for existing information and the application of theory.When studying a chemical reaction experimentally, some information such as the details of unstable intermediates and their rate constants are often inaccessible.To predict profiles of the intermediate steps, a plausible reaction mechanism must be proposed.To validate and verify the proposed mechanism and to estimate the unknown rate constants, simulations need to be done [6].For simple reactions involving one or two competitive or consecutive steps, estimation of the unknown rate constants is relatively easy process.When it involves multi-step mechanisms and a large number of rate constants to be estimated, it becomes a challenge.To estimate the unknown rate coefficients, the available data on similar type of reaction steps from the literature provides a guide line.Many of the radical chain reactions, branched chain reactions, chemical waves and Belousov-Zhabotinsky type reactions have complex reaction mechanisms [7].The validation of proposed mechanisms based on the simulations, has become an essential norm in assessing the complex chemical systems [8,9].While molecular modelling involves the prediction of structure of chemical compounds using computational chemical techniques, chemical modelling mostly deals with numerical simulation of chemical reactions with respect to time.Kinetic simulation is the imitation of the experimental behaviour with respect to time and involves the investigation and determination of the various reactant, product and intermediate species with respect to time and requires the proposed mechanisms for a reaction, the rate constants for each reaction step and the initial concentrations of the starting species [10].
Many of the commercial packages, such as NAG (Numerical Algorithms Group), ISML Inc, Numerical Recipes, provide subroutines in developing the kinetic simulation software [11].Mathematica, Matlab, etc. provide subroutines usable in developing appropriate software for simulations.The copyright, licensing and the high prices associated with that software restricts their flexible use.Hence freely accessible software with technically sound routines and mathematically compatible with chemical reactions, which can be duly custom-built, is great in demand.Most of the available software has been developed by programmers who are mostly mathematicians and have little or no exposure to all chemistry and depend on the information they received from the chemists.Software poses problems and difficulty to chemists who have limited or no experience in mathematics and programming.This program is developed with an objective of providing an user-friendly software which recognises and rectifies the user's errors, and to be able to simulate any complex chemical reaction mechanism with up to 40 reaction steps and can optimize the estimated rate constants for the chosen reaction steps within the specified or accepted range of values so that the simulated curves match with experimental curves.
We have earlier reported the programme "Simkine" and "Simkine2", and the theory and principles governing the development of that software [12,13].Initially "Simkine" was developed using Turbo Pascal for Windows 1.5 (TPW1.5)which was used to develop software for Windows 3.1 and Window 95.Simkine was further improved to give more flexibility and options and was rewritten using Borland Delphi Professional 5 (BD5) [14] and "Simkine2" which suits for Windows 2000, XP and Vista, and can operate at high speed, in addition to the other advantages of the programming language.In this article, we describe the software "Simkine3" in Delphi, which provides a prospect for optimisation of the rate coefficients for the chosen reactions to be estimated for the simulations.

THEORY
A chemical system pre-processor code is designed to translate a user specified system of chemical equations termed as mechanism in to a system of chemical kinetic differential equations termed as differential rate laws.The system of differential rate laws is translated into a system of integrated rate laws termed concentrations.In a mechanism each differential equation for a particular species depends on concentrations of other species and rate constants of each elementary step [15].A pre-processor code for a sequence of chemical reactions is described, which translates the user-specified system of chemical rate equations into a system of differential equations [11,16,17].The differential rate law in a mechanism can be represented as follows: symbolically, these differential equations can be represented as, ( ) The numerical solution of equation 2 is accomplished by using the semi-implicit midpoint rule and extrapolation as follows: the semi-implicit midpoint rule is used to approximate x(t+H) whenever x(t) is given, using m steps of size h = H/m.The result is denoted by x

(t + H, h). A polynomial extrapolation is applied to approximate x(t + H, h) by means of the data obtained from several values of h.
For practical implementation, if and let I denote the identity matrix in the following numerical scheme: For computation, the system becomes stiff as the time taken by some elementary steps differs due to the magnitudes of rate constants, i.e. they are consumed while the others are still in reaction.To handle the stiff problems, where computational efficiency is not crucial, the approach generally employed is semi-implicit Runge-Kutta method, which was used by Kaps and Rentrop in their program [18].The Predictor-Corrector methods, which are the Gear's backward differentiation methods and the generalizations of the Bulirsch-Stoer methods also called semi-implicit extrapolation methods (SIEM) [16].These are the mostly used methods for solving stiff problems in chemical systems.The SIEM methods gained popularity in recent past due to their increasing stability with increasing step size, and efficiency/speed with smaller error tolerance.The Simkine 2 program was developed using SIEM, which uses the modified midpoint rule.
With Simkine 3 the numerical solution was accomplished by using subroutine TLSoda developed by Sauro [19], a solver for first order differential equations (ODEs), because of its capability to solve stiff and non-stiff equations.The unit for optimisation of estimated constants contains the routines to perform non-linear regression to estimate the estimated rate constants so that the generated data matches the experimental data.Downhill Simplex Method (DSM) was used to perform non-linear regression.The code that implements DSM is taken from Numerical Recipes and the code for DSM subroutine which was originally written in Fortran was translated to Pascal without any modifications [11].
The subroutine STIBS, used in Simkine to solve the resulting differential equations (ODEs) is usual and typecast.When it comes to stiff problems, it takes a longer time to calculate the concentrations for all species that are present in the mechanism.As a result, it will hang up and may crash the system.In Simkine2 and Simkine3, instead of STIBS, the subroutine, TLSoda is employed for solving the generated ODEs.When considering the theory, the subroutine, TLSoda is also based on the functionality of STIBS, but it uses backward differentiation method to solve stiff problems [14].Instead of a MS Word file, an easy to read online help file written using 'HelpScribble' is provided as a separate file to assist the Simkine3 users [20].Furthermore, the necessity to store the data on the selected drive after each run is avoided, as it damages the concerned drive and may end up crashing the system.In 'Simkine3', the generated data is stored internally and can be saved later by option.
While simulating the kinetic profiles, user needs to view, how close the experimental and generated curves agree.Most of the available software during the simulations, plot the generated and experimental data separately and few have option to plot both for comparison at a later stage.Majority of the software allow the plots only after finishing the computation of the data.Overcoming these limitations, in order to compare the two curves, 'Simkine3' provides the plot of the generated and experimental data on the same set of axes simultaneously.Importantly, the present software provides the means to opt for the estimated rate coefficients to be optimized and again to set the ranges with in which those values to be optimized.These provisions allow the user to restrict the ranges with in which the estimated rate coefficients should be, based on the scientific experience and information available from the related studies.The optimization function will read experimental and generate data for the species that was varied while optimizing the estimated rate constants so that the generated data matches the experimental data.The estimated rate constants will be optimized using subroutine Amoeba which uses the Downhill Simplex Method [14].The experimental and generated data will be drawn on the same set of axes.Funk will be called to generate data for each curve, i.e. for each concentration for the varied species.

EXPERIMENTAL
The kinetic studies of the reaction between acidic bromate and methylene violet were conducted at (25.0 ± 0.1) o C with excess concentrations of acid and bromate and monitoring absorbance change of MV + (at 557 nm) using either Cary II-Varian UV-Visible spectrophotometer or the double mixing micro-volume stopped-flow (Hi-Tech Scientific SF-61DX2) spectrophotometer [22].Both the instruments were interfaced for data storage and have software to analyse the kinetic data.The reaction between methylene violet and aqueous bromine was investigated using the stopped-flow apparatus.The bromide ion concentration, during the progress of reaction, was monitored using the bromide ion selective (Orion) and a double junction reference electrode.

RESULTS AND DISCUSSION
The advantages of the 'Simkine3' software over the other software are (i) no programming experience is needed to use the program, (ii) it comes as the executable file ready to be run, therefore the user needs no contact with the source code, (iii) the user has to only supply the mechanism only as the text file using any editor that support text format, (iv) the reactant concentration conditions can be fed directly, (v) the number of the elementary steps and the species in the mechanism can be up to 40, (vi) all the data of the species available in the mechanism will be saved as a text file, which can be analysed using any spreadsheet that supports text format, (vii) the data and modified mechanism files can be stored or deleted by choice, (viii) the rate coefficients that are to be optimised and their ranges can be chosen and (ix) the user need not to be concerned about the fitting equation to model the data.A help file is provided to familiarise and to provide information on any queries.
To show the versatility of the 'Simkine3' program, a complex 19-step mechanism scheme, which was proposed for the oxidation of the heterocyclic phenazine dye, methylene voilet (MV + ) by bromate ion acidic conditions is chosen [26,27].
where MV + is methylene violet (3-amino-7-(diethyleamino)-5-phenyl phenazium ion and MP is 3,7-dioxo-5-phenyl phenazine [22].Academic interest in the methylene violet/acidic bromate reaction arises from its nonlinear kinetics, i.e. slow depletion rate at the initial stages followed by a very fast reaction after an induction time.During reaction, bromide ion, a reaction intermediate plays a dual role, both as the auto catalyst and as an inhibitor giving raise to the complex reaction mechanism.In the Belousov-Zhabotinskii type chemical systems, bromide ion, which switches between high and low concentration conditions acts as control intermediate [23,24].The change from high [Br -] to low occurs through the reaction of bromide with HBrO 2 and BrO 3 -.Regeneration of Br - depends upon the nature of the organic intermediates and their reaction rates with various oxybromo-and bromo-species in the system.The important equations representing the chemistry of key bromo and oxybromo intermediates in the reactions of bromate under acidic conditions [25][26][27] can be represented as follows:

R7
In the MV + and acidic bromate system, in the absence of initially added bromide, the induction time reflects the time required for the accumulation of HBrO 2 and the removal of trace inhibiting species, including bromide, present as impurities with bromate.The ratelimiting step in such system is direct reaction between MV + and BrO 3 -in presence of acid, which is followed by the formation of BrO 2 • radical [27,28].The reactions of organic intermediates with BrO 3 -, HBrO 2 , and HOBr are retarded by high [bromide] through its competitive reaction with the oxybromo species.High [bromide] inhibits the autocatalysis step (R6) by means of R2, depleting the bromous acid (HBrO 2 ) concentration.Thus, the accumulation of HBrO 2 is delayed, prolonging induction period for the rapid depletion of MV + .The fast depletion step is the result of the rapid reaction of MV + with Br 2 and HOBr.The increased concentration of HBrO 2 enhances its disproportion rate.At high [HBrO 2 ], with bromide concentration at critical level, bromide acts as auto catalyst (R3), resulting in an increase in [HOBr].Both Br 2 and HOBr could react rapidly with MV + and other organic species leading to the fast depletion of the organic substrate.Relative levels of HOBr and Br 2 are controlled by the concentrations of H + and Br -ions (R4).In addition, HBrO 2 also directly reacts with the organic intermediates [9,25].Based on the kinetic and product data, the following simplified steps may be proposed of the oxidation of MV + to MP.

Simulations
Using 'Simkine3', the characteristic features of the reaction were simulated.Literature values of rate coefficients are used for reactions, R1 to R7 [9,26,27] and for R19 [29].Experimental values of the rate constants are employed for R8 and R12 [22].The value of rate coefficient for R8 is rate-limiting and crucial.Estimated rate constants are used for the remainder of the reactions.Estimated rate coefficients for reactions involving intermediates with bromo and oxybromo species were kept high, and values were optimised such that the simulated curves agreed with the experimental curves [22].Table 1 summarizes all the elementary steps, the rate coefficients from the literature or experiments where available, with roughly assigned constants for the rate constants that have to be optimized.The estimated rate constants that need to be optimized are specified with a label '@" and the range within which the optimization should occur is laid down for each reaction.Table 2 shows an example of elementary reaction steps chosen for such optimization.Table 3 illustrates the roughly estimated rate constants and the ranges chosen with the acceptable minimum and maximum values for optimisation.The simulations can be repeated till the computed and experimental curves match well, to the user's approval.Finally, Table 4 shows the optimised rate constants for different reactions, which generate simulated curves in agreement with experimental curves.A help file is provided to familiarize and to provide information to any queries.Figure 1 illustrates the experimental curves and corresponding simulated curves for the variation of initial bromate showing that the simulated curves give good fits with all the six experimental curves simultaneously.The optimised values can be further refined, prior to conclusion of the simulation process.After the successful simulation, the concentration-time data for all the species for each set of reaction conditions generated can be saved as spread sheet.Figure 2 shows the experimental and simulated curves for [MV+] and simulated profiles of selected intermediates.The appearance of bromine only after complete depletion MV + agrees well with the reported fast reaction of bromine with substrate [27].Further, the effect of added bromide on the reaction profiles are also simulated (Figure 3).The fair agreement between the experimental and simulated curves in presence of varied amounts of added initial bromide further supports the validity of the proposed mechanism and near accuracy of the estimated rate coefficients.At low [Br -] 0 , simulations were sensitive to its change, the experimentally observed prolonged induction times with increase in added [Br -] could be fairly simulated.A good agreement in the induction times and profiles in the three sets of experimental and simulated curves (Figures 1 and 3) suggests that using the proposed reaction scheme and SIMKINE3, the optimization of the unknown rate coefficients could be achieved to acceptable harmony.Table 2. Rate coefficients to be optimized * .
Table 3. Rate coefficients to be optimized with an option to set the ranges * .
# The range of minimum and maximum values for the rate constant to be optimised can be set for each step.

CONCLUSIONS
This software offers an user-friendly, adaptable and tough version of the program to the chemistry fraternity, now branded as 'Simkine3' to compute any complicated multi-step chemical mechanisms and optimise an array of estimated rate constants with ease.