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

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.


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) Mathematics 2020, 8 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): 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].

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: 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 Mathematics 2020, 8, 1324 3 of 16 whose solutions are given by (Figure 1), whose solutions are given by (Figure 1), Note that in the first quadrant, due to the transformation of one species into another, the initial concentration of the species x (xstable + x0) tends to reduce to a local minimum (xstable), while species y tends to increase from its initial value (ystable) up to a local maximum value (ystable + y0). 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 (Ak1, k2, and k3), 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 (xstable, ystable).
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 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 0 ) 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 0 ). Consequently, rotation of the system will be counterclockwise, Figure 2. dt = dt = 0 (5) whose solutions are given by (Figure 1), Note that in the first quadrant, due to the transformation of one species into another, the initial concentration of the species x (xstable + x0) tends to reduce to a local minimum (xstable), while species y tends to increase from its initial value (ystable) up to a local maximum value (ystable + y0). 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 (Ak1, k2, and k3), 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 (xstable, ystable).
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 In the general case, it can be assumed that the interaction between species, which depend on the parameters of the problem (Ak 1 , k 2 , and k 3 ), 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):

1.
A value for x i above x stable (x i > x stable ) with y i = k 1 A/k 2 , 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.

2.
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.

3.
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.

4.
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: Final point: The dimensionless variables, which range in the interval [0, 1], are defined as follows: or, reorganizing, 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: 1 t o , −k 3 , k 2 x o , which provide the monomials (or dimensionless groups): On the other hand, from the dimensionless form of Equation (1), the emergent coefficients are 1 t o , k 1 A, and −k 2 y o , leading to the monomials 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 From these, applying the π-theorem, solutions for y o and t o are obtained, 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: Final point: Mathematics 2020, 8, x FOR PEER REVIEW 6 of 20 The dimensionless variables are defined in the form y = y y + k A k The dimensionless variables are defined in the form Proceeding as in case (i), Equation (1) provides the coefficients 1 t o , k 1 A, and −k 2 y o , and the monomials. while Equation (2) gives rise to the coefficients 1 t o , k 2 x o , −k 3 , and to the monomials.
The final monomials can be reorganized mathematically and written in the form providing the solutions Case (iii). The path is shown in Figure 4. Starting and final points are Starting point: Final point: Case (iii). The path is shown in Figure 4. Starting and final points are Starting point: Final point: These used as references provide the dimensionless variables These used as references provide the dimensionless variables Proceeding as in the former cases, the resulting monomials can be written in the form Thus, solutions for t o and x o are 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: Final point: while dimensionless variables are defined in the form while dimensionless variables are defined in the form Replacing these expressions in the governing equations, the resulting monomials are: Provide the solutions for t o and x o

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. 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: a.
Firstly, values are assigned to the parameters Ak 1 , k 2 , k 3 , x i , y i , and x o for cases (i) and (iii), and y o for cases (ii) and (iv). b.
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. Checked cases in Lotka-Volterra oscillator problem.

Parameters Variables
Case Scenario Case Scenario Case Scenario Case Scenario  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 (Figures 6-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 2 being halved. By contrast, in the second scenario, t o is doubled by the effect of halving k 3 , while y o does not change. Similar comments can be made for the rest of the simulations in the table. Figures 6-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 = −1 for cases (i) and (iii) and t o = 1 , −k 2 y o k 1 A 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: This fact can be checked immediately from This fact can be checked immediately from This fact can be checked immediately from Table 1 (Table 1). (a) Scenario X, (b) Scenario XI, (c) Scenario XII, and (d) phase plane of scenario X.  (Table 1). (a) Scenario X, (b) Scenario XI, (c) Scenario XII, and (d) phase plane of scenario X.

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.