Study of Lotka–Volterra Biological or Chemical Oscillator Problem Using the Normalization Technique: Prediction of Time and Concentrations

Abstract

The normalization of dimensionless groups that rule the system of nonlinear coupled ordinary differential equations defined by the Lotka–Volterra biological or chemical oscillator has been derived in this work by applying a normalized nondimensionalization protocol. The normalization procedure, which is quite accurate, does not require complex mathematical steps; however, a deep physical understanding of the problem is required to choose the appropriate references to define the dimensionless variables. From the dimensionless groups derived, the functional dependences of some unknowns of interest are established. Due to the coupled nature of the problem that induces temporal concentration rates of each species that are quite different at each point of the phase diagram, this diagram has been divided into four stretches corresponding to the four quadrants. For each stretch, the limit values (maximum or minimum) of the variables, as well as their duration, are expressed in terms of the dimensionless groups derived before. Finally, to check all the mentioned dependences, a numerical simulation has been carried out.

1. Introduction

The Belousov–Zhabotinsky reaction (BZ reaction) was initially discovered by Boris Belousov in the 1950s. The notable fact characterizing the BZ reaction, and it is extremely relevant because of this, is that the concentration of some species involved oscillates harmonically, showing typical structures of systems out of equilibrium. This reaction not only has great importance from the physical and chemical point of view, as well as mathematical, but it is also extremely suggestive from a biological point of view. The system of Lotka–Volterra differential equations belongs to these types of chemical or biological nonlinear oscillators, which has a special interest in the field of Systemic Thinking. In its simplest form, the Lotka–Volterra model is formed by two species that interact with one another, and its behavior reflects other models in different fields of science such as plasma physics [1], neural networks [2], ecology [3], chemistry [4,5,6], and biology [7,8].

By means of a simple mathematical manipulation [9,10,11], the normalized nondimensionalization reduces the governing equations to their dimensionless form, from which the dimensionless groups that rule the problem are derived. As is known, based on the π-theorem, these groups provide substantial information of the problem without solving it analytical or numerically [12,13]. This (Lotka–Volterra) kind of problem has a harmonic solution in which the total period can be split into four stretches whose durations are generally quite different, slowing down or speeding up the movement of the point in the phase-diagram cycle. These detailed unknowns, as well as the concentration limits that reach the variables at the end of each stretch, are functions of the aforementioned dimensionless groups (a direct application of the π-theorem).

Within the non-dimensionalization process, the suitable choice of references to make dimensionless and normalized the dependent and independent variables (which may be different for each one of the coupled equations) is perhaps the step that requires the most attention, as there are many possible parameters explicitly or implicitly included in the problem statement [11,12,13,14,15,16,17]. Some of these references are, in fact, unknowns whose order of magnitude is found once the dimensionless groups are deduced. In short, the suitable choice requires an exhaustive study up to a deep comprehension of the physical or chemical meaning of the terms involved both in governing equations and in boundary conditions.

Once the governing equations are made dimensionless, their terms are formed by two factors: One is a grouping of parameters; the other is a function of the dimensionless variables and their changes. Due to the normalization of the variables, the last factor is of an order of magnitude of unity in all the terms of the equation. In the balance between the pair of terms, the dimensionless groups are the direct ratios between factors, whose quotients are of the order of magnitude of unity. Through mathematical manipulation, we can sort the groups so that each unknown appears only in one of them. After that, if there are p dimensionless π_v-groups, each one containing a different unknown coefficient, and q dimensionless π_u-groups without unknowns, the solution for each p-group will be a determined arbitrary function (Ψ_i) of the q-groups (π-theorem): π_{v,i (1 ≤ i ≤ p)} = Ψ_i (π_u,1, π_u,2,… π_u,q).

If all the dimensionless groups (π) are of the order of magnitude of unity, the mentioned arbitrary function will also have this property. From the latter p-relations, the order of magnitude of each p-unknown is obtained.

The normalization procedure is illustrated with the Lotka–Volterra biological or chemical oscillator problem. To check the dimensionless numbers and the expressions that relate the unknown references to the rest of the parameters, a large number of cases have been carried out, numerically solving the problem using the software CODENET_15 [18], based on the Network Method methodology [19,20,21].

2. Lotka–Volterra Chemical Oscillator Problem

Using the nomenclature of Lotka [22] applied to the chemical model [23] or biological model [8], the governing equations are:

$$\frac{dx}{dt}=k_1{Ax-k}_{2}xy$$

(1)

$$\frac{dy}{dt}=k_2{xy-k}_{3}y$$

(2)

$$\frac{dA}{dt}=-k_{1}Ax$$

(3)

$$\frac{dB}{dt}=k_{3}y$$

(4)

In these equations, x and y are the main or intermediate species, while A and B are the reactant and product species, respectively. Assuming that the concentration of the A reactant, which generates the species x, Equation (1), is constant, the mechanism involved in these equations can be re-interpreted as a model with oscillating concentrations of intermediate products x and y. Focusing the study on these equations, Equations (1) and (2), at the equilibrium point, we can write

$$\frac{dx}{dt}=\frac{dy}{dt}=0$$

(5)

whose solutions are given by (Figure 1),

$$\begin{gathered} x_{stable}=\frac{k_3}{k_2} \\ {y}_{stable}=\frac{k_{1}A}{k_2} \end{gathered}$$

(6)

Note that in the first quadrant, due to the transformation of one species into another, the initial concentration of the species x (x_stable + x₀) tends to reduce to a local minimum (x_stable), while species y tends to increase from its initial value (y_stable) up to a local maximum value (y_stable + y₀). Consequently, rotation of the system will be counterclockwise, Figure 2.

In the general case, it can be assumed that the interaction between species, which depend on the parameters of the problem (Ak₁, k₂, and k₃), consists of a cycle of increase and decrease in the concentration of both, which will be periodically repeated. Temporal representations of these concentrations will be harmonic continuous waves, with different periods of rising and falling from its stable values (x_stable, y_stable).

Reducing Equations (1) and (2) to their dimensionless form, which is done following a normalized nondimensionalization process of the same, we obtain the independent dimensionless groups that, based on the π-theorem, govern the solutions of the problem. To do this, two references are needed for making dimensionless variables x and y. For simplicity, we choose initial values (x_i or y_i) located on the lines x = x_stable or y = y_stable. The possible alternatives for the choice of these references (x_i, y_i) are (subscripts i and f mean initial and final points, respectively, of the path segment in the phase diagram):

• A value for x_i above x_stable (x_i > x_stable) with y_i = k₁A/k₂, Equation (6). Then, the trajectory from the value (x_i, y_stable) to (x = x_stable,y_f) is studied. In this process, x_i is given and y_f is an unknown. The time for this trajectory is also unknown.
• A value for y_i above y_stable (y_i > y_stable) with x_i = x_stable. The trajectory goes from the value (x_stable, y_i) to the point (x_f,y = y_stable). In this process, y_i is given while x_f and the time required for the trajectory are unknowns.
• A value for x_i under x_stable (x_i < x_stable) with y_i = y_stable. The trajectory goes from the value (x_i, y_stable) to (x = x_stable,y_f). x_i is given while y_f and the time are unknows.
• A value for y_i under y_stable (y_i < y_stable) with x_i = x_stable. The trajectory goes from the value (x_stable, y_i) to the point (x_f,y = y_stable). y_i is given while x_f and the time are unknowns.

Different solutions for each of the previous cases will be provided because of no symmetry of the evolution of the species along the whole cycle. The problem is studied between two states of evolution contained in the lines x = x_stable and y = y_stable (one state in each line). The choice of an arbitrary starting point, where it is not contained in any of the lines x = x_stable or y = y_stable, is also possible. In this case, the starting point is (x_i, y_i) while the final, (x_f, y_f), is located at lines x = x_stable or y = y_stable, depending on the quadrant in which the starting point is located and the sense of rotation.

Let us deduce the dimensionless groups of the former alternatives:

Case (i). Locate at the first quadrant of the phase space (Figure 2). The starting point is the maximum value of the species x (distance x_o above the stable value x_stable). Write:

Starting point:

$$x_i=x_{max}=x_o+\frac{k_3}{k_2},y_i=y_{stable}=\frac{k_{1}A}{k_2}$$

(7)

Final point:

$$x_f=x_{stable}=\frac{k_3}{k_2},y_f=y_{max}=y_o+\frac{k_{1}A}{k_2}$$

(8)

The dimensionless variables, which range in the interval [0, 1], are defined as follows:

$$\begin{gathered} x^{\prime }=\frac{x-\frac{k_3}{k_2}}{x_{max}-\frac{k_3}{k_2}}=\frac{x-\frac{k_3}{k_2}}{x_o} \\ {y}^{\prime }=\frac{y-\frac{k_{1}A}{k_2}}{y_{max}-\frac{k_{1}A}{k_2}}=\frac{y-\frac{k_{1}A}{k_2}}{y_o} \end{gathered}$$

(9)

or, reorganizing,

$$\begin{gathered} x={x^{\prime }x}_o+\frac{k_3}{k_2} \\ y={y^{\prime }y}_o+\frac{k_{1}A}{k_2} \end{gathered}$$

(10)

Replacing these variables in the governing equations and defining t_o as the time elapsed between the initial and final points, the dimensionless forms of Equations (1) and (2) are deducted. On the one hand, from the dimensionless form of Equation (2), the following coefficients result: $\frac{1}{t_o}$, −k₃, k₂x_o, which provide the monomials (or dimensionless groups):

$$\begin{gathered} {\pi }_1=\frac{-1}{k_3{t}_o} \\ {\pi }_2=\frac{-k_2{x}_o}{k_3} \end{gathered}$$

(11)

On the other hand, from the dimensionless form of Equation (1), the emergent coefficients are $\frac{1}{t_o}$, k₁A, and −k₂y_o, leading to the monomials

$$\begin{gathered} {\pi }_3=\frac{1}{k_1{At}_o} \\ {\pi }_4=\frac{-k_2{y}_o}{k_{1}A} \end{gathered}$$

(12)

Now, manipulating groups (11) and (12) in order introduce the unknowns y_o and t_o in separated monomials, we can write them in the form

$$\begin{gathered} {\pi }_a=\frac{-k_{1}A}{k_3} \\ {\pi }_b=\frac{-k_2{y}_o}{k_{1}A} \\ {\pi }_c=\frac{-k_2{x}_o}{k_3} \\ {\pi }_d=\frac{1}{k_1{At}_o} \end{gathered}$$

(13)

From these, applying the π-theorem, solutions for y_o and t_o are obtained,

$$t_o=\frac{-1}{k_3}{\Psi }_1\left(\frac{-k_{1}A}{k_3},\frac{-k_2{x}_o}{k_3}\right)$$

(14)

$$y_o=\frac{-k_{1}A}{k_2}{\Psi }_2\left(\frac{-k_{1}A}{k_3},\frac{-k_2{x}_o}{k_3}\right)$$

(15)

Case (ii). The phase diagram is within the second quadrant, Figure 3. The starting point is the maximum value of the species y (y_o+y_stable). Write:

Starting point:

$$x_i=x_{stable}=\frac{k_3}{k_2},y_i=y_{max}=y_o+\frac{k_{1}A}{k_2}$$

(16)

Final point:

$$x_f=x_{min}=\frac{k_3}{k_2}-x_o,y_f=y_{stable}=\frac{k_{1}A}{k_2}$$

(17)

The dimensionless variables are defined in the form

$$\begin{gathered} x^{\prime }=\frac{x-x_{min}}{\frac{k_3}{k_2}-x_{min}}=\frac{x-x_{min}}{x_o} \\ \frac{y-\frac{k_{1}A}{k_2}}{y_{max}-\frac{k_{1}A}{k_2}}=\frac{y-\frac{k_{1}A}{k_2}}{y_o} \end{gathered}$$

(18)

or

$$\begin{gathered} x=x^{\prime }\left\{\frac{k_3}{k_2}-x_{min}\right\}+x_{min}={x^{\prime }x}_o+\left(\frac{k_3}{k_2}-x_o\right) \\ y={y^{\prime }y}_o+\frac{k_{1}A}{k_2} \end{gathered}$$

(19)

Proceeding as in case (i), Equation (1) provides the coefficients $\frac{1}{t_o}$, k₁A, and −k₂y_o, and the monomials.

$$\begin{gathered} {\pi }_1=\frac{1}{k_1{At}_o} \\ {\pi }_2=\frac{-k_2{y}_o}{k_{1}A} \end{gathered}$$

(20)

while Equation (2) gives rise to the coefficients $\frac{1}{t_o}$, k₂x_o, −k₃, and to the monomials.

$$\begin{gathered} {\pi }_3=\frac{-1}{k_3{t}_o} \\ {\pi }_4=\frac{-k_2{x}_o}{k_3} \end{gathered}$$

(21)

The final monomials can be reorganized mathematically and written in the form

$$\begin{gathered} {\pi }_a=\frac{-k_{1}A}{k_3} \\ {\pi }_b=\frac{-k_2{y}_o}{k_{1}A} \\ {\pi }_c=\frac{-k_2{x}_o}{k_3} \\ {\pi }_d=\frac{1}{k_1{At}_o} \end{gathered}$$

(22)

providing the solutions

$$t_o=\frac{1}{k_{1}A}{\Psi }_3\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$$

(23)

$$x_o=\frac{-k_3}{k_2}{\Psi }_4\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$$

(24)

Case (iii). The path is shown in Figure 4. Starting and final points are

Starting point:

$$x_i=x_{min}=\frac{k_3}{k_2}-x_o,y_i=y_{stable}=\frac{k_{1}A}{k_2}$$

(25)

Final point:

$$x_f=x_{stable}=\frac{k_3}{k_2},y_f=y_{min}=\frac{k_{1}A}{k_2}-y_o$$

(26)

These used as references provide the dimensionless variables

$$\begin{gathered} x^{\prime }=\frac{\frac{k_3}{k_2}-x}{\frac{k_3}{k_2}-x_{min}}=\frac{\frac{k_3}{k_2}-x}{x_o} \\ {y}^{\prime }=\frac{\frac{k_{1}A}{k_2}-y}{\frac{k_{1}A}{k_2}-y_{min}}=\frac{\frac{k_{1}A}{k_2}-y}{y_o} \end{gathered}$$

(27)

or

$$\begin{gathered} x=\frac{k_3}{k_2}-{x^{\prime }x}_o \\ y=\frac{k_{1}A}{k_2}-{y^{\prime }y}_o \end{gathered}$$

(28)

Proceeding as in the former cases, the resulting monomials can be written in the form

$$\begin{gathered} {\pi }_a=\frac{-k_{1}A}{k_3} \\ {\pi }_b=\frac{-k_2{y}_o}{k_{1}A} \\ {\pi }_c=\frac{-k_2{x}_o}{k_3} \\ {\pi }_d=\frac{1}{k_1{At}_o} \end{gathered}$$

(29)

Thus, solutions for t_o and x_o are

$$t_o=\frac{1}{k_{1}A}{\Psi }_5\left(\frac{-k_{1}A}{k_3},\frac{-k_2{x}_o}{k_3}\right)$$

(30)

$$y_o=\frac{-k_{1}A}{k_2}{\Psi }_6\left(\frac{-k_{1}A}{k_3},\frac{-k_2{x}_o}{k_3}\right)$$

(31)

Case (iv). The phase space locates at the fourth quadrant, Figure 5. The starting point is the minimum value of the species y (y_stable − y_o). The references are now:

Starting point:

$$x_i=x_{stable}=\frac{k_3}{k_2},y_i=y_{min}=\frac{k_{1}A}{k_2}-y_o$$

(32)

Final point:

$$x_f=x_{max}=\frac{k_3}{k_2}+x_o,y_f=y_{stable}=\frac{k_{1}A}{k_2}$$

(33)

while dimensionless variables are defined in the form

$$\begin{gathered} x^{\prime }=\frac{x_{max}-x}{x_{max}-\frac{k_3}{k_2}}=\frac{x_{max}-x}{x_o} \\ {y}^{\prime }=\frac{\frac{k_{1}A}{k_2}-y}{\frac{k_{1}A}{k_2}-y_{min}}=\frac{\frac{k_{1}A}{k_2}-y}{y_o} \end{gathered}$$

(34)

or

$$\begin{gathered} x=x_{max}-x^{\prime }\left\{x_{max}-\frac{k_3}{k_2}\right\}=\left(\frac{k_3}{k_2}+x_o\right)-{x^{\prime }x}_o \\ y=\frac{k_{1}A}{k_2}-y^{\prime }{y}_o \end{gathered}$$

(35)

Replacing these expressions in the governing equations, the resulting monomials are:

$$\begin{gathered} {\pi }_a=\frac{-k_{1}A}{k_3} \\ {\pi }_b=\frac{-k_2{y}_o}{k_{1}A} \\ {\pi }_c=\frac{-k_2{x}_o}{k_3} \\ {\pi }_d=\frac{1}{k_1{At}_o} \end{gathered}$$

(36)

Provide the solutions for t_o and x_o

$$t_o=\frac{1}{k_{1}A}{\Psi }_7\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$$

(37)

$$x_o=\frac{-k_3}{k_2}{\Psi }_8\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$$

(38)

3. Verification of Results

In order to check the accuracy and reliability of the expressions deduced for unknowns in the four paths of the former section, we numerically simulate the set of twelve scenarios shown in

Table 1. Checked cases in Lotka–Volterra oscillator problem.

				Parameters							Variables
Case	Scenario	xst = k3/k2	yst = k1A/k2	k1	k2	k3	A	xi	xo	yi	yo	to
I	I (Basis)	1	1	1	1	1	1	2	1	1	0.998	1.196
	II (2·to)	1	1	0.5	0.5	0.5	1	2	1	1	0.998	2.296
	III (2·yo)	2	2	1	0.5	1	1	4	2	2	1.996	1.196
Case	Scenario	xst = k3/k2	yst = k1A/k2	k1	k2	k3	A	xi	yo	yi	xo	to
Ii	IV (Basis)	1	1	1	1	1	1	1	0.5	1.5	0.374	1.596
	V (0.5·to)	1	1	2	2	2	1	1	0.5	1.5	0.374	0.796
	VI (2·xo)	2	2	1	0.5	1	1	2	1	3	0.748	1.596
Case	Scenario	xst = k3/k2	yst = k1A/k2	k1	k2	k3	A	xi	xo	yi	yo	to
Iii	VII (Basis)	1	1	1	1	1	1	0.5	0.5	1	0.500	2.096
	VIII (2·to)	1	1	0.5	0.5	0.5	1	0.5	0.5	1	0.500	4.196
	IX (2·yo)	2	2	1	0.5	1	1	1	1	2	1	2.096
Case	Scenario	xst = k3/k2	yst = k1A/k2	k1	k2	k3	A	xi	yo	yi	xo	to
Iv	X (Basis)	1	1	1	1	1	1	1	0.5	0.5	0.756	1.596
	XI (0.5·to)	1	1	2	2	2	1	1	0.5	0.5	0.756	0.796
	XII (2·xo)	2	2	1	0.5	1	1	2	1	1	1.512	1.596

. These scenarios are grouped in four cases, with three scenarios in each case. The first scenario of each case is chosen as a basis in order to be compared to the other two simulations. The procedure may be summarized in two steps:

• Firstly, values are assigned to the parameters Ak₁, k₂, k₃, x_i, y_i, and x_o for cases (i) and (iii), and y_o for cases (ii) and (iv).
• Secondly, these values are suitably changed in such a way that the dimensionless groups involved in each case retain the same value. Therefore, the arguments of the functions from which the unknowns t_o and y_o (or x_o) depend do not change, and neither do the functions themselves.

Table 1 shows the values of the parameters of each scenario. Solutions for x_o and t_o in cases (ii) and (iv), respectively, and y_o and t_o in cases (i) and (iii), respectively, read from the simulations (Figure 6, Figure 7, Figure 8 and Figure 9), are included in the last columns of the table. Note that, for example, in case (i), y_o duplicates in the third scenario (in comparison to the first scenario) due to the parameter k₂ being halved. By contrast, in the second scenario, t_o is doubled by the effect of halving k₃, while y_o does not change. Similar comments can be made for the rest of the simulations in the table. Figure 6, Figure 7, Figure 8 and Figure 9 show the harmonic evolution of species and allow t_o to be read for each scenario. Phase diagrams are also shown in the figures.

It is clear that the dependences between unknowns and parameters (or monomials), Equations (14) and (15) for case (i), (23) and (24) for case (ii), (30) and (31) for case (iii), and (37) and (38) for case (iv), leads to a different time t_o for each one of the quadrant, a solution that clearly emerges from the simulations. Note, however, that the solutions for t_o can also be made in the same form for the four cases by simple mathematical manipulation. In effect, the solutions $t_o=\frac{-1}{k_3}{\Psi }_1\left(\frac{-k_{1}A}{k_3},\frac{-k_2{x}_o}{k_3}\right)$ for cases (i) and (iii) and $t_o=\frac{1}{k_{1}A}{\Psi }_3\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$ for cases (ii) and (iv) are interchangeable. This means that the total period of a complete cycle of the problem (T_o), the sum of the four times of each quadrant, is also dependent on the same monomials:

$$T_o=\frac{1}{k_{1}A}{\Psi }_9\left(\frac{-k_{1}A}{k_3},\frac{-k_2{y}_o}{k_{1}A}\right)$$

(39)

This fact can be checked immediately from Table 1 and Figure 6, Figure 7, Figure 8 and Figure 9.

4. Final Comments and Conclusions

The normalized nondimensionalization of a system of two coupled ordinary differential equations (that represent problems of interactions between species type Lotka–Volterra) has been a relatively quick and reliable procedure that allows one to obtain the dimensionless groups that rule the solution of these problems. From these groups and applying the π-theorem, the order of magnitude of some unknowns of interest can also be determined. To transform the governing equations into their dimensionless form, we only need simple mathematical manipulations; however, to choose the appropriate reference values that make the independent and dependent variables of equations dimensionless, an in-depth understanding of the physical or chemical processes involved has been required.

A general conclusion can be observed regarding the large variety of possible references for the independent and dependent variables. The references for a specific variable should not necessarily be the same in each equation of the system nor in the analysis of the solution of each quadrant of the phase diagram. In the Lotka–Volterra system under study, we have several references for concentrations as there are initial values, as well as the stable point concentrations. In addition, normalization of the dimensionless variables, which range in the interval [0, 1], imposes new restrictions for the choice of reference values. Furthermore, due to the nonlinearity of the problem, the characteristic time has been defined for each portion of the quadrant of the diagram phase in order to obtain the dependence of these partial times on the dimensionless groups of the problem. Even if the times depend on the same groups, the unknown functions of the arguments (monomials) are not necessarily the same, so the times are of the same order of magnitude but not of the same value. In addition, the period corresponding to a whole cycle (which obviously depends on the same dimensionless groups) has not been obtained directly but as a sum of the times related to the four quadrant trajectories. Overall, it is shown that the procedure used in this work gives rise to very valuable information about the solutions of a complex phenomenon, just with a little mathematical effort.

In view of the results, the normalized nondimensionalization procedure applied in this work opens two perspectives: (a) The construction of universal abacuses that represent the unknowns of interest as a function of the dimensionless groups, on the one hand (a result derived from the π- theorem); for this end, a large number of numerical solutions is required. (b) The application of the procedure to systems with more than two variables, on the other hand. However, as the number of parameters increases significantly with the number of variables, the complexity of the study increases greatly due to the existence of different periods in each phase diagram.