Theoretical Foundation for the Discrete Dynamics of Physicochemical Systems:
Chaos, Self-Organization, Time and Space in Complex Systems
V. Gontar
International Group for Scientific and Technological Chaos Studies, Ben-Gurion University of the Negev, P.O.B. 653, Beer-Sheva 84105, Israel

 

Abstract

A new theoretical foundation for the discrete dynamics of physicochemical systems is presented. Based on the analogy between the p-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

    Since Lorenz's discovery of chaotic solutions of differntial equations [1], H. Haken's analysis of self-organized systems [2], I. Prigogine's discussion about chaos and order, time and space in complex non-equilibrium systems [3] , and L. Chua's chaotic electrical circuirts [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 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, synchronize the output of electronic circuits, control oscillations in chemical reactions, stabilize the erratic 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 tiny, intricate patterns from which models of virtually anything, from snowflakes to mountains, can be created. The insights from this exercise can be so exciting that Wall Street firms are funding extensive research into fractals and its sister field, chaos theory, to see if they can predict stock market behavior." (Business Week, p. 169, October 12, 1992).
    "Motile cells of Escherichia coli aggregate to form stable patterns of remarkable regularity when grown from a single point on certain substrates. Central to this self-organization is chemotaxis, the motion of bacteria along gradients of a chemical attractant that the cells themselves excrete. . . . Formation of spatial patterns from a mass of initially identical cells is one of the central problems of developmental biology." (Nature, v. 376, no. 6535, 6 July 1995).
    "The discipline of chaos has created a universal paradigm, a scientific parlance, and a mathematical tool for grappling with nonlinear phenomena. In every field of the applied sciences (astronomy, atmospheric sciences, biology, chemistry, economics, geophysics, life and medical sciences, physics, social sciences, zoology, etc.) and engineering (aerospace, chemical, electronic, civil, computer, information, mechanical, software, telecommunication, etc.), the local and global manifestations of Chaos and Bifurcation have burst forth in an unprecedented universality, linking scientists heretofore unfamiliar with one another's fields, and offering an opportunity to reshape our grasp of reality." (Aims and Scope, International Journal of Bifurcation and Chaos. World Scientific).
    It is obvious that we are faced with a new dynamic paradigm demanding explanation and much better understanding. How this new paradigm of chaos exist alongside with the classical dynamics, based on the use of differential equations for which infinitesimal changes of initial conditions and parameters result in infinitesimal changes in solutions? Chaotic solutions, which in opposition to the conventional solutions of differential equations, are extremely sensitive to infinitesimal changes - bringing into question the main postulate of the correct use of differential equations. Remembering the fact that chaotic solutions can be obtained only by computer, we must ask whether they reflect the dynamic laws of nature (written in the form of differential equations before chaotic regimes were proposed), or whether they are solely the result of an extreme sensitivity to numerical procedures and computer errors [5]. How can the continuous time and space of differential equations reflect catastrophic changes in behavior in chaotic systems, and how can we combine the determinism and predictability of differential equations and the unpredictability of chaotic solutions? What is the new role of computers (which have become an essential part of our solutions of differential equations through the accuracy of numerical procedures) when chaos appears and changes in accuracy can strongly affect the results of calculation? It is obvious that computing errors never can be made "as small as necessary"—one of the main postulates of all existing numerical methods and the ideology of using computers—if chaotic solutions are expected. The old contradiction between continuality and discreteness is thus once again clamoring for an answer. The traditional chain of differential equations - numerical solution has thus been broken by the simple fact that chaotic solutions depend not only on differential equations but strongly on numerical algorithms and computer errors. A new link between classical mathematical models and their solutions should be put into the chain. This link should include the accuracy of numerical methods, computer accuracy, etc.
    There are a minimum of two ways of solving these problems. The first one is to make a careful analysis of the numerical accuracy of integration of differential equations (which is very hard, and not always possible, to do). The second one is to seek other methods (calculus) different from differential equations and free of the problems connected with their use. The question whether systems of ordinary or partial differential equations are the single mathematical language for describing all kinds of dynamics derived from first principles of physics is not new. Can iterational 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 calculus—calculus of iterations? If we succeed with this program, we not only overcome many contradictions that result from the use of differential equations, but we also take a step towards the better understanding of the meaning of discrete time & space, self-organization & complexity. Of course, any activity 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 results. This is the program that we intend to undertake, 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 p-theorem of the theory of dimensionality, which bring us to the formulation of a new extreme 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 application to the description of different complex multicomponent systems. The principle of a minimum of free energy for closed multicomponent systems and following mass action law of thermodynamics provides us with the necessary mathematical models in a form algebraic equations for the case of equilibrium:
(1)
and differential equations of the kinetic mass action law for the non-equilibrium case:
{ – 
(2)
where Kl - the equilibrium constant for the l reaction, and  - rate constants for the l reaction, li- matrix of stochiomeric coefficients, D - diffusion constants.
    Despite the absence of a variational principle for the rate equations (2), they are considered as fundamental equations that have verified by numerous experiments in chemistry, physics, biology, etc. In the case that equation (2) is applied to simulations of chaotic processes, 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 nli. For all types of constituents Ai (atoms, molecules, radicals, ions, clusters, cells, etc.), the mechanism of interaction can be written in the following form:

 
nli Ai = 0    i = 1, 2, . . . , N
   l = 1, 2, . . . , N – M 
(3)
    When nli is given as an initial hypothesis, all the necessary equations for simulation can be automatically written according to the equations (1) or equations (2). By adding the equations of the law of mass conservation:
aij Xi = bj         j = 1,2, . . . , M
(4)
we will complete the system of algebraic equations in the case of equilibrium (1) and the system 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 nli and aij the following relationship should be satisfied:
(5)
This relation can be used in order to calculate matrix aij from nli, or matrix nli, when matrix aij 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 dimension expressed with the aid of the so-called main unitsmeters, seconds, kilograms, etc. Variables such as velocity, energy, pressure, etc. have well-defined dimensions, which are obviously dependent on the choice of the main units. To overcome the dependence of arbitrary choice of the main units, we need, and recommend, to construct invariants (i.e., dimensionless variables). Invariants are constructed by using the p-theorem of the theory of dimensionality and are more conveniently expressed by the general form of matrix operation [9].
    Suppose we have N variables qi with the dimensions defined by use of matrix Qij, where j = 1, 2,. . . , and M - number of main units; for example, q1 = s - distance, q2 = t - time, q3 = v - velocity, q4 = a - acceleration, variables defining the kinematics of an arbitrary mechanical system with main units: u1- sec (s), and u2 - meter (m). In this case, matrix Q can be written in the form:
 
m
 s
s
|
1  0 |
t
|
0  1 |
v
|
1 -1 |
a
|
1 -2 |
(6)
Now to define the invariants according to the p-theorem, we need to present matrix Qij in the following form:
 
|
R
|
Q =
|
– – – –
|
|
P
|
 
(7)
where R is nondegenerate matrix of dimensions M ´ M, and P is a matrix of dimensions (N-M) ´ M. Then, the determinant matrix zli for L=N-M invariants pl take the following form:
zli = |-PR-1 I |
(8)
where I is the unit matrix (N-M) ´ (N-M).
For our example we will get matrix:
(9)
and two invariants:
p1 = s-1 t1 v1 a0 = vt/s
 
p2 = s-1 t2 v0 a1 = at2 /s.
(10)
The general expression for N-M invariants pl is given by the p-theorem:
 
(10a)
Now let us transport this formalism to the stoichiometry of chemical reactions.
Again, to make things clear, let us use a simple example of a chemical system with four constituents X1 = A, X2 = B, X3 = AB, X4 = AB2 and as the components of the system, let us choose A and B. In this case, the molecular matrix aij takes the following form:
 
A
B
 
A
|
1
0
|
B
|
0
1
|
AB
|
1
1
|
AB2
|
1
2
|
 
(11)
and the chemical "invariants" are easy to define in the same way used in equation (8). As a result we will get the following expressions for the two "invariants":
(12)
    Now, the analogy between the equations of the thermodynamic mass action law (1) and the p-theorem (10a) is obvious: the matrixes aij & Qij and nli & zli are mathematically equivalent, as well as equations (1) and (10a). Of course, there is a vast difference in the meaning of the dimensionless invariants pl and the "chemical invariants" Pl, which have a dimension in the sense of the theory of dimensionality, but nonetheless we intend to use some similarities and the analogies.
    As is known, equations (1) were derived in classical thermodynamics from the principle of a minimum of free energy (or maximum entropy) for a mixture of ideal gases. This assumption serves as a strong constraint for the use these equations for a condensed systems. If we view this from the position of the theory of dimensionality and suppose that the dimensionless character of the "chemical invariants" reflects the independence of the thermodynamic constants Kl (mathematical analog to pl) on the concentrations of M components of the system, we can apply equations (1) to any systems with chemical reactions, without constraints, 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 situation 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 zli (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 physiochemical systems, we need to add a new "time component" in the molecular matrix aij, when going to dynamic simulations from thermodynamic equilibrium. The details of such an approach can be found in ref. [10]. If we do this for a closed system, we will get, after simple operations of the matrixes, the time-dependent functions Pl(tq) = Kl exp[-Wl/tq] on the right side of equations (1), where q = 1,2. . . :
= Kl exp[-Wl/tq]
(13)
which together with equations (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 asymptotic behavior of chemical systems [10]. The equations (4,13) ha 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) without using of differential equations (2) for any mechanisms of chemical reactions given by matrix nli. 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 mechanisms of chemical reactions, and good qualitative 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 approximation to the solutions of differential equations.
    Continuing to the mathematical simulations of open catalytic systems, we intend to include in our molecular matrix aij a new component reflecting time delay interaction, which after simple transformations [10], will bring us to a system of difference equations. In this case, Pl 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 t0 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 coordinate r through dependence Xi(rÄ) in our matrix aij, and the final function for Pl can be written in the form:
(14)
where ali, bli are empirical parameters, characterizing the intensity of feedbacks in time - and the intensity influences of the space-distributed concentrations of the neighbors Xi(rÄ) - on the lth 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 nli.
Using the fact that equations (1) express the minimization 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 considered space r(r1, r2, r3) the function:
F(tq,r) =     {lnXi(tq) – fi – 1} i = 1,2, . . . ,N
(15)
reaches its minimum in the concentration space Xi, subject to the to the law of mass conservation (4).
Here:
 

where Pl is defined by equation (14).

Results of Numerical 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 mechanism of chemical transformations:

(16)
where A = X1, B = X2, C = X3 and X3 affects reaction A®B and X1 affects reaction B®C.
    We constrain ourselves with the case, in which only the previous moment of time is included in the consideration s=1. According to the mechanism (15), we have the following stoichiometric and molecular matrixes:
 
 
|
-1
1 0 |
|
1
|
nli
= | |      
aij
=
|
1
|
|
-1
0
1
|
|
1
|
 
 
and equations (13,4):
 
 
 
X2(tq,r) 
--------
=
K1{[a13X3(tq–1) + b11X1(tq–1,rÅ) + b12X2(tq–1,rÅ) + b13X3(tq–1,rÅ)]/tc}
X1(tq,r)
X3(tq,r)
--------
=
K2{[a11X1(tq–1) + b21X1(tq–1,rÅ) + b22X2(tq–1,rÅ) + b23X3(tq–1,rÅ)]/tc}
X1(tq,r)
 
 
X1(tq,r) + X2 (tq,r) + X3 (tq,r) = b

(17)

    Concerning spatiotemporal dynamics, we constrain ourselves in this publication to the plane r = r(r1, r2) represented by a square lattice with Â*Â elements with coordinates r(r1, r2). In each cell, we calculate Xi(tq,r) according to equations (17), which are reflect the influence of the concentrations of the neighbors Xi(tq-1,rÅ) on the reaction rate in each cell with coordinates r(r1,r2). The initial conditions are: Xi(tq,r) = b, X2(tq,r) = X3(tq,r) = 0. The boundary conditions reflect the fact of the absence of neighbors for all cells with coordinates |r|= Â.
    Our intention with this mathematical model (17) was to demonstrate the ability of the proposed theoretical approach to generate different types of complex chemical oscillations behavior and to simulate pattern formation dynamics.
In Picture 1, we present different time series obtained from equation (17).
Pictures 2 and 3 present a process of creation of a patterns that are growing and moving in space. Picture 4 demonstrates the creation of two patterns and their interaction. Picture 5 demonstrates the evolution of chemical spiral waves. In Picture 6 we present some symmetrical patterns. All patterns were obtained in a square lattice with 160´160 cells, except picture 4: 70x70 cells. The color of each cell reflects the value of the concentration of Xi(tq,r). Calculations were performed on a PC-486 computer.
 
Discussion and Conclusions
    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 equations from the computational point of view are obvious (the dynamics of pattern formation presented here takes about 10 min 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 problem of defining the meaning of so-called "discrete time and space" appearing when iterational maps or difference equations are used. The questions arise when we try to make a link between the continuous time and space of differential equations and discrete nature of our numerical calculations. As a precondition for sampling, we suppose the existence of a "small-as-necessary" time interval Dt. This concept fails when chaotic solutions can be obtained from differential 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 extremely 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, equations (1) and (13) demonstrate a new time tq. Trajectories Xi(tq) are defined by the solutions of the system of non-linear algebraic equations (1,4,13). Despite of the close similarity of these solutions to the solutions of differential equations, the meaning of time tq and that of continuous time t in (2) is different from the calculation point of view, and tq can be considered as discrete time by virtue of the way that this time appears. The continuality of time t in differential equations is connected to the existence of lim dX/dt and dt ® 0, in contrast to tq, which runs in arbitrary way from t0 in any order we need.
    The next step in the proposed theory is the description of chemical reactions with feedbacks, when the parameters of the initial mathematical model get the dependence of the concentration Xi calculated at the previous "moment of time tq-s". In this case, the algebraic equations transform into a system of difference equations, and for an open system our discrete time tq becomes a constant or in another 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 derivatives dX/dt) or even in the sense of the discrete time tq in equations (13,1,4). Astronomic continuous time t commonly used for interpretation 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 another words, that the system is "self-organizing in our physical time and space," and only the character of the interactions is responsible for any time intervals appearing in the system's dynamics. The smallest time interval ?t, 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 integrating 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 l 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 Picture 2, at the early stages of distribution of the concentration we see that concentrations are not organized in any specific form. But later, due solely to chemical transformations going in each cell of the square lattice we can see some process of concentration organization 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 obtain from the new basic equations is the direct concentration distribution of particular constituent 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 background [12].
    We should also relate to the stochastic process of mathematical simulation. Because chaotic regimes 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 existence of other theoretical and mathematical foundations to chemical reaction dynamics than the differential equations of the kinetic mass action law.
    We are perched at the beginning of the development 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 calculus of the infinitesimal for differential equations. The results already coming out from proposed approach prove the tremendous potential of the proposed mathematical model that can be effectively used in real-time control systems, neural networks, and in signal processing as practically unlimited source of different types of oscillatory 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 application in economics, medicine, computer science and many other theoretical and practical areas for which extensive computing time is required, without losing the physicochemical sense and theoretical background of mathematical models.
 
Pictures:
Pic.1
 

 
1
2
3
4
5
6
 
Pic.2 Dynamics of pattern formation, obtained by new basic equations for the chemical reactions mechanism: A-->B-->C. 1- 50th iteration, 2-100th, etc.

 

1
2
3
4
5
6
 
Pic.3 Dynamics of pattern formation for the same mechanism of chemical reactions, but different parameters of the mathematical model.


 
 Pic.4. Dynamics of patterns formation and their interaction.

 
1
2
3
4
5
6
 
 Pic.5. Evolution of spiral chemical waves , obtained by new basic equations for B-Z reaction: 1- 300th iteration, 2- 400th, ..., 6- 800th iteration.

 
  
Pic.6. Symmetrical patterns obtained by new basic equations.

Acknowledgment
I would like to express my gratitude to the Ben-Gurion University of the Negev, the Lady Davis Trust, and the Doron Foundation for Education and Welfare, for their support.
 
References
1. E.N. Lorenz, "Deterministic non periodic flow", J. Am. Sci. 1963, 20, 130-141
2. H. Haken, Synergetic: An Introduction, Springer Ser. Synergetics, Vol.1, 1978.
3. I. Prigogine, “From Being to Becoming: Time and Complexity in the Physical Sciences”, W.H. Freeman, San-Francisco, 1980.
4. L. O. Chua, C. A. Desoer, and E. S. Kuh, "Linear and Nonlinear Circuirts" New York, NY: McGraw-Hill, 1987.
5. J.H.E. Cartwright and O. Piro, “The Dynamics of Runge-Kutta Methods”, Int. Journal of Bifurcation and Chaos, V. 2, N. 3, 1992, pp. 427-450.
6. V. Gontar, “Theory of Dimensionality and Chemical Reactions Kinetics. A New Model of Chemical Kinetics”, Vest. Mosc. University, Chem. V. 17, N. 1, 1976, pp. 108-110.
7. V. Gontar, “New Principle and Mathematical Model for Chemical Reactions Kinetics”, R. Journal Phys. Chem. V. 55, N. 9, 1981, pp. 2301-2304.
8. V. Gontar, I. Ilin. “New Mathematical Model of Physicochemical Dynamics”. Contrib. Plasma Phys. 31, N. 6, 1991, pp. 681-690
9. L. Brand, “p-Theorem of the Theory of Dimensionality”, Arch. Rational Mech. Anal., V. 1, 1957, pp.35-43.
10 V. Gontar, “New Theoretical Approach for Physicochemical Reactions Dynamics with Chaotic Behavior”, In: “Chaos in Chemistry and Biochemistry”, World Scientific, London, 1993, pp. 225-247.
11. I. R. Epstein, K. Showalter, "Non-linear Chemical Dynamics: Oscillations, Patterns, and Chaos", J. Phys. Chem. 1996, 100, pp. 13132-13147.
12. T. Toffoli, N. Margolus: Cellular Automata Mashines - A New Environment For Modelling, MIT Press Cambridge/Mass, London/ Engl. !987.