Theoretical Foundation for the Discrete Dynamics of Physicochemical Systems : Chaos , Self-Organization , Time and Space in Complex Systems

A new theoretical foundation for the discrete dynamics of physicochemical systems is presented. Based on the analogy between the π-theorem of the theory of dimensionality, the second law of thermodynamics and the stoichiometry of complex physicochemical reactions, basic dynamic equations and an extreme principle were formulated. The meaning of discrete time and space in the proposed equations is discussed. Some results of numerical calculations are presented to demonstrate the potential of the proposed approach to the mathematical simulation of spatiotemporal physicochemical reaction dynamics.


INTRODUCTION
Lorenz's discovery of chaotic solutions of differential equations [1], Haken's analysis of self-organized systems [2], Prigogine's discussion about chaos and order, time and space in com- plex non-equilibrium systems [3], and Chua's chaotic electrical circuits [4], there have been so many important contributions to the field of complex non-linear chaotic system dynamics that it is impossible to do them all justice.There is no other field in science that could compete with such an interest and activity.More than two hundred books with the word chaos in the title 31 have been published, seven new journals devoted to chaos studies have been launched in the last few years... "What is good in chaos?... Chaos has already been applied to increase the power of lasers, syn- chronize the output of electronic circuits, control oscillations in chemical reactions, stabilize the er- ratic beat of unhealthy animal hearts and encode electronic messages for secure communications.We anticipate that in the near future engineers will no longer shun chaos but will embrace it."(Scientific American, p. 78, August 1993).
"Fractals are an effort to simulate nature's complexity.Those computer representations are V. GONTAR tiny, intricate patterns from which models of vir- differential equations, are extremely sensitive to tually anything, from snowflakes to mountains, infinitesimal changes-bringing into question the can be created.The insights from this exercise main postulate of the correct use of differential can be so exciting that Wall Street firms are equations.Remembering the fact that chaotic so- funding extensive research into fractals and its lutions can be obtained only by computer, we sister field, chaos theory, to see if they can pre- must ask whether they reflect the dynamic laws dict stock market behavior."(Business Week,  of nature (written in the form of differential p. 169, October 12, 1992).
equations before chaotic regimes were proposed), "Motile cells of Escherichia coli aggregate to or whether they are solely the result of an ex- form stable patterns of remarkable regularity treme sensitivity to numerical procedures and when grown from a single point on certain sub- computer errors [5].How can the continuous strates.Central to .thisself-organization is chemo- time and space of differential equations reflect taxis, the motion of bacteria along gradients of a catastrophic changes in behavior in chaotic sys- chemical attractant that the cells themselves ex- terns, and how can we combine the determinism crete... Formation of spatial patterns from a mass and predictability of differential equations and of initially identical cells is one of the central the unpredictability of chaotic solutions?What is problems of developmental biology."(Nature, the new role of computers (which have become v. 376, no. 6535, 6July 1995).
an essential part of our solutions of differential "The discipline of chaos has created a universal equations through the accuracy of numerical pro- paradigm, a scientific parlance, and a mathemati- cedures) when chaos appears and changes in cal tool for grappling with nonlinear phenomena, accuracy can strongly affect the results of cal- In every field of the applied sciences (astronomy, culation?It is obvious that computing errors atmospheric sciences, biology, chemistry, eco- never can be made "as small as necessary"one nomics, geophysics, life and medical sciences, of the main postulates of all existing numerical physics, social sciences, zoology, etc.) and engi- methods and the ideology of using computers- neering (aerospace, chemical, electronic, civil, if chaotic solutions are expected.The old contra- computer, information, mechanical, software, tele- diction between continuality and discreteness is communication, etc.), the local and global mani- thus once again clamoring for an answer.The festations of Chaos and Bifurcation have burst traditional chain of differential equations-numerforth in an unprecedented universality, linking ical solution has thus been broken by the simple scientists heretofore unfamiliar with one an- fact that chaotic solutions depend not only on other's fields, and offering an opportunity to re- differential equations but strongly on numerical shape our grasp of reality."(Aims and Scope, algorithms and computer errors.A new link be- International Journal of Bifurcation and Chaos.tween classical mathematical models and their World Scientific).solutions should be put into the chain.This link It is obvious that we are faced with a new should include the accuracy of numerical meth- dynamic paradigm demanding explanation and ods, computer accuracy, etc. much better understanding.How this new para- There are a minimum of two ways of solving digm of chaos exists alongside with the classical these problems.The first one is to make a careful dynamics, based on the use of differential equaanalysis of the numerical accuracy of integration tions for which infinitesimal changes of initial of differential equations (which is very hard, conditions and parameters result in infinitesimal and not always possible, to do).The second one changes in solutions?Chaotic solutions, which is to seek other methods (calculus) different from in opposition to the conventional solutions of differential equations and free of the problems connected with their use.The question whether systems of ordinary or partial differential equa- tions are the single mathematical language for describing all kinds of dynamics derived from first principles of physics is not new.Can itera- tional maps or systems of difference equations be used instead of or in parallel with differential equations?What type of postulates do we need to accept in order to put a sound physicochemical meaning to this new calculuscalculus of iterations?If we succeed with this program, we not only overcome many contradictions that re- sult from the use of differential equations, but we also take a step towards the better understanding of the meaning of discrete time and space, self- organization and complexity.Of course, any ac- tivity in finding the required principles should be based on the traditional procedure of verification: equations derived from the new principles should not only describe existing experimental results, should not contradict solutions of differential equation for non-chaotic regimes, and should predict a wide.range of future experimental re- suits.This is the program that we intend to un- dertake, and with this publication, which presents one of the possible approaches, we open the discussion.
Starting from classical chemical thermodynamics for closed systems, we have found an analogy between the second law of thermodynamics and the 7r-theorem of the theory of dimensionality, which bring us to the formulation of a new ex- treme principle for chemical reaction kinetics [6,7] and dynamics [8].

BACKGROUND
The formalism of using the notations of chemical reactions is a very common and effective applica- tion to the description of different complex multi- component systems.The principle of a minimum of free energy for closed multicomponent systems and following mass action law of thermody-namics provides us with the necessary mathema- tical models in a form of algebraic equations for the case of equilibrium: and differential equations of the kinetic mass ac- tion law for the non-equilibrium case: ( where KIthe equilibrium constant for the /th reaction, and k-, krate constants for the lth reaction, uli matrix of stoichiomeric coefficients, Di-difflusion constants.Despite the absence of a variational principle for the rate equations (2), they are considered as fundamental equations that have been verified by numerous experiments in chemistry, physics, biology.In the case that Eq. ( 2) is applied to simulations of chaotic pro- cesses, the same questions that were put forward in the Introduction must be raised.
To initiate the procedure of constructing a new mathematical model, we need to start by defining the mechanism of chemical reactions, expressed by matrix of stoichiometrical coefficients lli.For all types of constituents Ai (atoms, molecules, radicals, ions, clusters, cells, etc.), the mechanism of interaction can be written in the following form: When lli is given as an initial hypothesis, all the necessary equations for simulation can be auto- matically written according to Eqs. (1) or Eqs.(2).By adding the equations of the law of mass V. GONTAR conservation: N aiXi bj, j 1,2,..., M, (4) i=1 we will complete the system of algebraic equa- tions in the case of equilibrium (1) and the sys- tem of differential equations (2) in the case of chemical kinetics.Here aij-"molecular" matrix, defining the number of components of type "j" in the ith constituent; and Xi-concentration of ith constituent of the system.According to the stoichiometrical rules of chemical transformations, for any matrixes uli and a the following relationship should be satisfied: This relation can be used in order to calculate matrix a O. from l/li or matrix lli when matrix a is given.Now let us turn to the foundation of the theory of dimensionality.According to the claims of this theory, all the variables used to identify the physicochemical systems should have a di- mension expressed with the aid of the so-called main units-meters, seconds, kilograms, etc. Variables such as velocity, energy, pressure, etc. have well-defined dimensions, which are ob- viously dependent on the choice of the main units.To overcome the dependence of an arbi- trary choice of the main units, we need, and recommend, to construct invariants (i.e., dimen- sionless variables).Invariants are constructed by using the 7r-theorem of the theory of dimension- ality and are more conveniently expressed by the general form of matrix operation [9].
Suppose we have N variables qi with the di- mensions defined by use of matrix Qo, where j--1,2,..., and Mnumber of main units; for example, ql =s-distance, q2--t-time, q3 u- velocity, q4 a-acceleration, variables defining the kinematics of an arbitrary mechanical system with main units: Ul sec (s), and U2 meter (rn).In this case, matrix Q can be written in the form Now to define the invariants according to the r-theorem, we need to present matrix Q0 in the following form: R (7 where R is non-degenerate matrix of dimensions M xM, and P is a matrix of dimensions (N-M) x M.Then, the determinant matrix li for L=N-M invariants r take the following form: )i=l_pR-1 ii, where I is the unit matrix (N-M) x (N-M).
For our example we will get matrix -1 0 -1 2 0 (9) and two invariants: 7r --s -1 l/1 aut/s, 7r2 S-t2 l/0 al atZ/s. (10) The general expression for N-M invariants rz is given by the 7r-theorem: N 7r, H q,i. Now let us transport this formalism to the stoi- chiometry of chemical reactions.Again, to make things clear, let us use a simple example of a chemical system with four constituents X A, X2 B, X3 =AB, X4 =AB2 and as the compo- nents of the system, let us choose A and B. In this case, the molecular matrix a O. takes the fol- lowing form: and the chemical "invariants" are easy to define in the same way used in Eq. ( 8).As a result we will get the following expressions for the two "in- variants": Now, the analogy between the equations of the thermodynamic mass action law (1) and the theorem (10a) is obvious: the matrixes aij and Q0 and uti and li are mathematically equivalent, as well as Eqs.( 1) and (10a).Of course, there is a vast difference in the meaning of the dimension- less invariants 7rt and the "chemical invariants" 1-It, which have a dimension in the sense of the theory of dimensionality, but nonetheless we in- tend to use some similarities and the analogies.
As is known, Eqs.(1) were derived in classical thermodynamics from the principle of a mini- mum of free energy (or maximum entropy) for a mixture of ideal gases.This assumption serves as a strong constraint for the use of these equations for a condensed system.If we view this from the position of the theory of dimensionality and sup- pose that the dimensionless character of the "chemical invariants" reflects the independence of the thermodynamic constants Kt (mathematical ana- log to 7rt) on the concentrations of M compo- nents of the system, we can apply Eqs.(1) to any systems with chemical reactions, without con- straints, based on ideal gases assumptions.If the initial system of chemical reactions does not give an adequate description, according to the theory of dimensionality, we need to change our initial hypothesis about the mechanism and add some more components and constituents, or change the type of chemical reaction.
Another recommendation that can be taken from the theory of dimensionality is the claim for an invariant to be constant (not dependent on time, initial concentrations, or concentrations of constituents) for similar systems.If a new situa- tion arises (for example, from simulation of a kinematic system, we are going to analyze the dynamics), we need to add a new component and variables in the matrix li (kg as a main unit in this particular case, and force, energy, as new variables etc.).Following this rule for moving to the description of the dynamics of the physico- chemical systems, we need to add a new "time component" in the molecular matrix ao., when going to dynamic simulations from thermo- dynamic equilibrium.The details of such an ap- proach can be found in Ref. [10].If we do this for a closed system, we will get, after simple op- erations of the matrixes, the time-dependent functions IIt(tq)= Kt exp[Wt/tq] on the right side of Eqs.(1), where q 1,2...: N H Xi'i(tq K exp[-W/tq], (13) 1-1 which together with Eqs.(4) represent a complete system of N non-linear algebraic equations for N variables Xi(tq).The reasons for our choice of an exponential function can be found in the asymp- totic behavior of chemical systems [10].Eqs.(4), (13) have a unit solution for all Xi(tq) > 0. This is the way in which we can construct a system of algebraic equations for mathematical simulation of the dynamics of N concentrations Xi(tq) with- out using differential equations (2) for any mechanisms of chemical reactions given by V. GONTAR matrix l/li.Numerous numerical calculations have been made to compare the solutions that have come out of the ordinary differential equations (without diffusion) of the kinetic mass action law and the proposed equations for the same me- chanisms of chemical reactions, and good quali- tative agreement was obtained [6].System (13), ( 4) is far more simple from a computational point of view (no problems with the stiffness of differential equations) and can serve as good ap- proximation to the solutions of differential equations.
Continuing to the mathematical simulations of open catalytic systems, we intend to include in our molecular matrix a/ a new component re- flecting time delay interaction, which after simple transformations [10], will bring us to a system of difference equations.In this case, IIl becomes a function of all concentrations Xi, calculated at previous moments of time tq-s, s--1,2,...For an open system, when tq changes from to to tc- a constant (characteristic time for any specific system, when the steady state is reached), system (13),(4) transforms from algebraic equations to a set of non-linear difference equations.In the same conceptional way, we insert space coordi- nate r through dependence Xi(r(R)) in our matrix aij, and the final function for IIl can be written in the form IIl(X(tq),r) Klexp{-N Z Oli i(tq-s) i=1 + Z tli Xi(r@) 'c (14) i=1 where Oli /li are empirical parameters, character- izing the intensity of feedbacks in time and the intensity influences of the space-distributed con- centrations of the neighbors Xi(r(R))on the /th chemical reaction.
We have now constructed a full system of basic equations for simulations of all types of complex chemical reaction dynamics, starting with the initial hypotheses of the mechanism given by matrix b'li.
Using the fact that Eqs.(1) express the mini- mization of free energy, we would like to generalize this principle and formulate a new extreme principle for chemical reaction dynamics: Reactions in multicomponent physicochemical systems proceed in such a way that at each instant of time tq (q--1,2...) and at every point of consid- ered space r(rl,r2, r3) the function: N F(tq, r)-Z Xi(tq, r){ln Xi(tq, r)-f.-1},i=l i= 1, 2,...,N, (15) reaches its minimum in the concentration space Xi, subject to the to the law of mass conservation (4).

CALCULATIONS
Some results of numerical calculations can now be presented to demonstrate possible applications of our approach to the mathematical simulation of the spatial dynamics of pattern formation and time series.Let us consider the following me- chanism of chemical transformations: T where A X1, B X2, C X3 and X3 affects re- action A --B and X1 affects reaction B-+ C.
The initial conditions are: Xi(tq, r) b, Xz(tq, r) X3(tq, r) --0.The boundary conditions reflect the fact of the absence of neighbors for all cells with coordinates Ir[ 9t.
Our intention with this mathematical model (17) was to demonstrate the ability of the pro- posed theoretical approach to generate different types of complex chemical oscillations behavior and to simulate pattern formation dynamics.In Fig. 1, we present different time series obtained from Eq. (17).Figures 2 and 3 present a process of creation of patterns that are growing and moving in space.Figure 4 demonstrates the crea- tion of two patterns and their interaction.Figure 5 demonstrates the evolution of chemical spiral waves.In Fig. 6 we present some symmetrical patterns.All patterns were obtained in a square lattice with 160 x 160 cells, except Fig. 4:70 x 70 cells.The color of each cell reflects the value of We have presented a brief description of a new theoretical approach to physicochemical reaction dynamics.We consider this approach to be the basis for constructing a solid theoretical foundation for the application of a system of algebraic and difference equations, instead of, or in parallel to, the commonly used system of differential equations.The advantages of difference equa- tions from the computational point of view are obvious (the dynamics of pattern formation pre- sented here takes about 10min for 300 iterations on a PC-486).Now let us clarify the meaning of time and space that we face in our approach and equations.We think that there is a general pro- blem of defining the meaning of the so-called "discrete time and space" appearing when itera- tional maps or difference equations are used.The questions arise when we try to make a link be- tween the continuous time and space of differen- tial equations and discrete nature of our numerical calculations.As a precondition, for sampling, we suppose the existence of a "smallas-necessary" time interval At.This concept fails when chaotic solutions can be obtained from dif- ferential equations.When problems of computer accuracy and error diffusion problems are also considered, an analysis of the real meaning and quality of the obtained solutions becomes extre- -mely complicated.
Our approach from the very beginning was directed towards constructing algebraic equations and then difference equations which are, by definition, free of the above-mentioned problem.For a closed system, Eq. (13) demonstrate a new time tq.Trajectories Xi(tq) are defined by the so- lutions of the system of non-linear algebraic equations (4), (13).Despite the close similarity of these solutions to the solutions of differential equations, the meaning of time tq and that of

FIGURE
Different oscillatory curves, obtained using new basic equations ( 4), ( 13), (14).,, Symmetrical patterns obtained using new basic equations.continuous time in (2) is different from the calculation point of view, and tv can be consid- ered as discrete time by virtue of the way that this time appears.The continuality of time in differential equations is connected to the exis- tence of lira dX/dt and dt --+ O, in contrast to tq, which runs in arbitrary way from to in any order we need.The next step in the proposed theory is the description of chemical reactions with feedbacks, V. GONTAR when the parameters of the initial mathematical model get the dependence of the concentration Xi calculated at the previous "moment of time tq_".In this case, the algebraic equations trans- form into a system of difference equations, and for an open system our discrete time tq becomes a constant or in other words disappears as a dynamic variable.Difference equations, being dynamic or evolutionary equations by their very origin, have no time in the sense of time of the systems of differential equations (2) (through de- rivatives dX/dt) or even in the sense of the dis- crete time tq in Eqs. ( 13), ( 1), (4).Astronomic continuous time commonly used for interpreta- tion of all types of dynamic behavior has disappeared from our difference equations, being transformed in the set of internal times resulting from our calculations.The initial scale of time is not important for difference equations and can be easily attached to the astronomic time in our experiments by normalization of parameters of the model when trajectory X is calculated.This means, in other words, that the system is "self- organizing in our physical time and space", and only the character of the interactions is responsi- ble for any time intervals appearing in the sys- tem's dynamics.The smallest time interval At, as well as all the other intervals, can be calculated by spectrum analysis of trajectories X.The fact that the smallest interval no longer depends on the limitations of the numerical procedure of in- tegrating differential equations is very important for practical computer calculations.
The situation with simulation of spatial dynamics by our system of difference equations is also different from use of partial differential equations (2).Spatial coordinates are not included in the basic equations, and are present only through dependence of parameters IIt from concentrations Xi(tq, r).The absence of derivatives in our basic equations frees us from the necessity of existence of continuous curves, as happened in the case of solutions of partial differential equations (2).
If we look at Fig. 2, at the early stages of distribution of the concentration we see that con- centrations are not organized in any specific form.But later, due solely to chemical transfor- mations going in each cell of the square lattice we can see some process of concentration organi- zation in a form that reminds us of spiral and ring curves.In fact, these are not curves in the sense of continuous curves that can be obtained from partial differential equations.What we ob- tain from the new basic equations is the direct concentration distribution of particular constitu- ent X in space according to our extreme principle and initial hypothesis of the chemical reaction mechanism.This principle is self-sufficient for the simulation of different types of spatiotemporal behavior of chemical waves and other pattern formations, which was previously the only partial differential equations [11] or "cellular automata" method with very little physicochemical back- ground [12].
We should also relate to the stochastic process of mathematical simulation.Because chaotic re- gimes are present in the proposed basic difference equations, there is no special need to add a "random number generator", as should be done with the differential equations.The random number generator and chaotic regimes are the same in their mathematical nature.
In summary, we have demonstrated the exis- tence of other theoretical and mathematical foun- dations to chemical reaction dynamics than the differential equations of the kinetic mass action law.
We are perched at the beginning of the devel- opment of this new paradigm for discrete dynamics and see as our goal for future studies the development of a calculus of difference equations in the same manner as it was done with the cal- culus of the infinitesimal for differential equations.The results already coming out from proposed approach prove the tremendous poten- tial of the proposed mathematical model that can be effectively used in real-time control systems, neural networks, and in signal processing as practi- cally unlimited source of different types of oscilla- tory and spatiotemporal behavior simulations.These approaches can be effectively used for mathematical simulations of living systems.We are sure that the application of the concept of proposed discrete dynamics and the new basic equations will find applications in economics, medicine, computer science and many other theo- retical and practical areas for which extensive computing time is required, without losing the physicochemical sense and theoretical back- ground of mathematical models.

FIGURE 2
FIGURE 2 Dynamics of pattern formation, obtained using new basic equations for the chemical reaction mechanism: A B--, C. represents the 50th iteration; 2 the 100th; etc.

FIGURE 3 FIGURE 4
FIGURE 3 Dynamics of pattern formation for the same chemical reaction mcchanism as in Fig.2, but using different para- meters in the mathematical model.

FIGURE 5
FIGURE5 Evolution of spiral chemical waves, obtained using new basic equations for B-Z reaction: corresponds to the 300th iteration; 2 to the 400th; 6 to the 800th.