Discrete Dynamics for Mathematical Simulation of Living 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 approach and mathematical model for simulation of the dynamic behavior of complex systems are presented. Chemical reaction dynamics are used as a basis for the mathematical modelling of the dynamics of living systems. The meaning of discrete time and space appearing in the new equations are discussed. Some numerical results of a new mathematical model for simulation of the oscillatory and spatial dynamic behavior of physicochemical reactions are given. The proposed theoretical and mathematical tools are considered as the beginning of the development of a "calculus of iterations" and difference equations in addition to the "calculus of infinitesimal" and differential equations for mathematical modelling of nonlinear dynamics.

 

Introduction

The age-old question about the possibilities of the mathematical simulation of living systems continues to arouse great scientific interest [1-3].

What are the main differences, if any, between the dynamics of living and nonliving matter from the mathematical point of view? How do the new achievements of modern mathematics, physics and biology affect this problem? How we can define living systems in mathematical terms? What are the consequences of the ideas coming out of chaotic systems analyses.

It is well known that—as opposed to nonliving systems—living systems can reproduce themselves and can collect and exchange information.

Chemical reaction kinetics, as expressed by the differential equations (DE) of the kinetic mass action law (KMAL), is the most widely used mathematical model for simulation of complex molecules, organisms, or population growth dynamics. The greater the volume of new experimental results about oscillations, self-organization procceses and other complexities in physicochemical systems [4], the more complex and non-linear the KMAL equations become, and the greater the computer time required to get adequate results[5]. The neccesity to include feedbacks in KMAL drastically increases the "stiffness" of DE and causes numerical and understanding problems, especially when chaotic regimes occur. Chaotic regimes—being by definition extremely sensitive to infinitely small errors and therefore to numerical procedures—often obfuscate the value of conclusions obtained from computers [6]. Do these conclusions really reflect our initial hypotheses about the mechanisms of interactions or are we faced with computer-added errors? If we change only the accuracy and method of numerical calculations, will we change our conclusions about the real physicochemical system? Can we achieve, in principle, the accuracy required for simulation of complex nonlinear systems with chaotic regimes? All these problems force us to look for other mathematical tools and theories that would be much simpler from a calculation point of view. At the same time, we are looking for new mathemathical approaches and theories encompassing general physicochemical principles and laws. In the consideration of modern non-linear sciences, it appears that we have a deficit of general ideas that could be applied to explain the questions discussed above and that current thinking has to be supplemented by new principles. The main goal of this paper is thus to demonstrate one possible new method for the mathemathical simulation of complex living systems.

 

Background

From our point of view, the dynamics of complex chemical reactions is the closest advanced language for our purpose—the mathematical simulation of living systems. Reaction kinetics of chemical transformations suppose that as a result of interactions between the constituents of any specific system—atoms, molecules, clusters or ions (or likewise, men, women, cities, etc.)—a new constituent or a new set of conditions appear in the system under consideration. Chemical reactions, by their very nature, include the fact of a creation, which is one of the most important features of living systems.

Let us denote all the constituents of the considered system by Ai. Then, we can formally represent a complex multicomponent physicochemical transformation in the form:

 

              i = 1,2, . . . , N 
l = 1,2, . . . , L
 
(1)
 
where nil - matrix of stoichiometric coefficients reflecting the number constituents of the type i in the lth chemical reaction, N - total number of constituents, and L - total number of chemical reactions, L = N – M.
For multicomponent systems, we should use the law of mass conservation:
 
              i=1,2, . . . , M
 
(2)
 
where aij - "molecular" matrix defining the quantity of the jth component in the ith constituent, bj - total amount of the jth component, and Xi - concentration of the ith constituent.
Components and constituents having the same nature can both belong to the list given by Ai, except that the total quantity of the jth component should be constant value in a closed system, according to eqn. (2). The total number of components M are given by the investigator as a hypothesis and may be changed during the process of mathematical model identification, as we do with the types and numbers of constituents as well. The more precise our experimental data and initial information about types and numbers of constituents, the more constituents and components can be identified and included in the mathematical model for simulation.
As was mentioned above, living systems should be differentiated from nonliving systems by another important feature, such as the possibility of information exchange. Of course, when we say information exchange, we should constrain ourselves by some concrete definition. In our particular case, by "information exchange," we intend to include in our calculations some special mathematical relationships between the concentrations Xi and the parameters of the mathematical model. In KMAL, the parameters are the rate constants, which by definition should not be dependent on the concentrations of the constituents. Nevertheless, for catalyzed chemical systems, there are some methods to define the relationships between concentrations and constants, which we will not discuss here, but instead we will present a new model, where establishing of this type of relationship looks natural. Thus, by "information exchange", we imply special types of relationships between the parameters of the mathematical model and the variables Xi. This means that we use the term "information exchange" simply to underline the difference between mass transformations given by relations (1) and some other transformations of informations, taking place in a system, by changing parameters of the mathematical model in a special way .
For this reason, we suggest the following way of formalization of the information exchange processes together with mass exchanges. To illustrate this idea let us take the simple bimolecular chemical reaction:
 
A + B « AB

(3)
 

Written in traditional way this notation reflects the mass exchange between the constituents A, B, and AB. To include the "information exchange" into mathematical formalism, we have to include the direct dependence of the concentrations of the constituents Xa, Xb, and Xab on the parameters responsible for the chemical reaction dynamics.

Formally, the influence of the constituents can be represented in the following , for example, form:

 

(4)
 

According to eqn. (4), as the initial hypotheses, constituent A accelerates the rate of the reaction, and constituent AB, in opposition, decreases it. In general the suitable form of chemical reaction mechanism of mass (what are the amounts and types of components and constituents?) and "information" (what are and how the constituents affecting the reaction) exchange should be defined from a statistical analysis of the experimental data. At this point we are ready to raise a question about the type of eqfor chemical reaction kinetics. At present the basis of some physicochemical foundation, the only kind of equations used, are the set of ordinary or partial differential equations of KMAL. However, recently obtained results in non-linear sciences have given us the confidence to claim that there are other mathematical and physical ideas reasonable enough to displace the differential equations from monopolizing the description of all kinds of dynamic systems and to clarify and make more reasonable the formalization of the "information exchange." The starting point for this approach was the intentions of computer scientists to decrease time of numerical calculation of differential equations. Continuous time and space of differential equations and a discrete type of computing is the driving force to look for the new theoretical ways for the mathematical simulation of complex nonlinear systems. The practice of developing new generations of computer systems cannot be considered as the only way to increase our computational facilities. Only together with the development of new physicochemical and mathematical approaches can we expect positive results in the analysis and simulations of complex system.

The theory we are suggesting for simulations of chemical reaction dynamics is based on an analogy to the p-theorem of the theory of dimensionality [7] and stoichiometry of chemical reactions [9]: 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 the considered space r(x,y,z) the function

 

              i = 1,2,. . . ,N
 
(5)
 
reaches its minimum in the concentration space Xi, subject to the to the law of mass conservation "(2).
Here

where pl are the "chemical invariants" [8] and the functions dependent on time tq, space coordinates r(x,y,z), temperature T, and pressure P. For closed chemical systems the following time dependence was proposed for pl

 

p l = K l exp [ -Wl/tq] l = 1, 2,....., L = N - M

(6)
 

where Kl are equilibrium constants and Wl are empirical parameters characterizing the rate of the lth reaction. Its obvious that when tq tends to infinity, the function F transfers to the regular free energy of multicomponent chemical systems.

Using this principle, we obtain a set of L nonlinear algebraic equations for simulation of chemical reaction dynamics [10]. The equivalent task to finding the minimum of F under the constrains (2) is to solve the following system of algebraic equations:

 

(7)

together with eqn. (2). Solving the system (2,7) for each moment in time tq , we obtain the dependences for N variables X i (tq ) as functions of time tq and parameters Wl, Kl , and bj and mechanism of chemical reactions, given by matrix nli.

Calculations made to compare solutions of the differential equations of the KMAL and suggested equations (2,6,7) showed good qualitative agreement for all possible mechanisms of chemical reactions in closed systems. This provided us with the backing to offer algebraic equations for simulation of the dynamics in closed chemical systems as a good, but much more simpler from the calculation point of view, approximation of the KMAL.

Returning to simulations of living systems, we need to expand our approach and formalize "information exchange" to take into consideration the case when the system is open and to include space coordinates in our equations. For these reasons ,we claim that our parameters Wl should have a dependence on the concentrations Xi calculated at previous instants of time tq-s, s = 1,2, . . . and on the concentration of the "neighbors" Xi(r). In our case, we want to expand the meaning of the word "neighbors" to include in the consideration all the constituents of the system which may be situated far from the particular point where the chemical reaction is taking place at that moment, but nevertheless affects the reaction. Not going into detail about different types of dependencies for Wl, we propose the following form for the dependence W on Xi(tq-s) and r:

Wl (Xi (tq–s), (r)) =(r)

(8)

where ali, bli are empirical parameters. In an open system tq changes from 0 to the constant time tc , when the steady state is reached and also should be considered as an empirical parameter, defined from the experimental data.

Now, we can complete construction of the discrete dynamics equations for the mathematical simulation of chemical reaction dynamics . They take the form of a system of nonlinear difference equations and are fully defined by a matrix of stoichiometric coefficients nli /the matrix aij is connected to nli by the matrix equation: nli aij = 0 and therefore may be easily determined [8]/. This means that for any complex system, whose behavior can be represented by the mechanism of physicochemical reactions transformation with mass and "information" exchange given by the matrixes ali, bli we can construct set of difference equations (2,7) with (6,8), solve them, and get as a result the time dependences and space distributions for all constituents Xi(tq,r).

These dynamic equations derived from the new extreme principle formulated here may be considered a theoretical approach that can contribute to our understanding of the meaning of time and space. In opposition to the continuous time and space in differential equations, we are faced here with "discrete" time and space. In addition to the tremendous advantages for computer calculations of the system of difference equations, we want to emphasize the special meaning of time that appears, when difference equations are used. Astronomic, continuous time disappears from our equations, and a set of "internal times", as a characteristic of the temporal behavior of the complex multicomponent system, emerges from our consideration. We are now not constrained by our watches as the only one scale of time, but a set of "internal times" will be obtained as a result of our calculations. (We may scale our obtained "times" to the astronomic time when needed). The same can be said about the space dynamics of complex systems—the geometrical forms that result from the calculations of our equations can and should be subsequently used in the regulary way, after ordinary scaling to real experimental space. Space patterns will reflect the mass and "information exchange" according to the particular mechanism of the reactions taking place and will represent forms according to our extreme principle (5).

Results of Numerical Calculations

To demonstrate the possibilities of the new discrete dynamics for mathematical simulations of complex system behavior we present some results of calculations made for the simple mechanism of two chemical reactions:

 

 
 

Xa  
Xb
Xc 
| -1 1 0 | |1|  Xa
nli  = | |
aij
=
|1|   Xb Xa + Xb + Xc=b
| 0 -1 1 |     |1|  Xc 
 
 
 
 

where concentration Xa affects reaction B ® C, and Xc influences reaction A ->B. We constraining our simulations by using equations (2) and by considering the case in which only the previous moment of time is included: s=1, and only close neighbors are included.

In picture 1 we present different types of temporal behavior for this mechanism of the chemical reaction and different values of the parameters b, K1,K2, ali, and bli. In Picture 2 the dynamics of pattern formation are presented. The colors correspond to the value of concentration X1 at particular points on the drawings. Each drawing consists of 160 x 160 cells and the initial concentrations Xi in each cell are equal: X1 = b, X2 = 0, X3 = 0. Picture 2a presents the 350th step of calculations, Picture 2b the 400th, Picture 2c the 500th, 650th, 750th, and 800t. The entire calculation took 4 min of conputer time on a Pentium PC, 166 MHZ.

 

Conclusions

In the proposed theoretical approach, we consider a mathematical tool for numerical modeling of dynamics of complex multicomponent systems.

As may be seen from some results presented in this paper, the new dynamic model has practically unlimited facilities to simulate any type of oscillations, including chaotic . Oscillations and waves play a very special important role in practically all disciplines of science. Chemical oscillations are a comparatively new subject for study and we think that new mathematical tools for their simulations are needed. One of the possible ways is presented here. We found good agreement between our theoretically calculated types of oscillatory curves and experimental results obtained in B-Z reaction, CO2-oxidation and ECG analyses.

We hope that a new calculus—"calculus of iterations", should be further developed and used in a same way as we do with calculus of infinitesimal and differential equations. Being much simpler from a calculation point of view (in comparison with differential equations) the proposed mathematical model can be used when real-time calculations are needed.

Our model is deterministic—and does not use any probabilistic approaches—for simulating the dynamics of chemical reactions in space and time. The set of difference equations, by their very nature, provide chaotic regimes, which play the same role as the introduction of the random number generator in differential equations, and therefore suitable for simulations of a system's complex stochastic behavior.

We understand that we are standing at the very beginning of the new discrete dynamics of complex and living systems, but even the preliminary results we have obtained by the proposed basic algebraic equations and subsequently difference equations with clear physicochemical meaning have given us the confidence to believe that we have entered a field with great potential.

 

  Pictures:

 
Pic. 1

 
 
a
 
b
 
c
 
d
 
e
 
f
 
Pic.2.  Spatial dynamics of chemical reactions A==B==C, calculated by new dynamics equations.
 
  
 
  References
1.I. Prigogine, “From Being to Becoming: Time and Complexity in the Physical Sciences”, W.H. Freeman, San-Francisco, 1980.
2. H. Haken and H.P. Koepchen, “Rhythms in Physiological Systems”, Proceedings of the International Symposium at Schloss ELMAU, Bavaria, October 22-25, 1990. Berlin, Springer, 1991.
3. M. Mitchell Waldrop, “Complexity: The Emerging Science at the Edge of Order and Chaos” Penguin Books, Harmondsworth, England, 1994.
4. E. O. Budrene and H. C. Berg, Dynamics of formation of symmetrical patterns by chemotactic bacteria, Nature, 376, 6535, 1995.
5. F. Argoul, A. Arneodo, P. Richetti, and J.C. Roux, “From Quasiperiodicity to Chaos in the Belousov-Zhabotinskii Reaction”, Journal Phys. Chem. V.76, N. 6, 1987.
6. J.H.E. Cartwright and O. Piro, “The Dynamics of Runge-Kutta Methods”, Int. Journal of Bifurcation and Chaos, V.2, N.3, 1992.
7. L. Brand, “p-theorem of the Theory of Dimentionality”, Arch. Rational Mech. Anal., V. 1, p.35, 1957.
8. V. Gontar, “Theory of Dimensionality and Chemical Reaction Kinetics. A New Model of Chemical Kinetics”, Vest. Mosc. University, Chem. V. 17, N. 1, 1976.
9. V. Gontar, “New Principle and Mathematical Model for Chemical Reactions Kinetics”, R. Journal Phys. Chem. V. 55, N. 9, 1981.
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.