A Nonlinear Structure of a Chemical Reaction Model and Numerical Modeling with the New Aspect of Existence and Uniqueness

: In this article, a nonlinear autocatalytic chemical reaction glycolysis model with the appearance of advection and diffusion is proposed. The occurrence and unicity of the solutions in Banach spaces are investigated. The solutions to these types of models are obtained by the optimization of the closed and convex subsets of the function space. Explicit estimates of the solutions for the admissible auxiliary data are formulated. An elegant numerical scheme is designed for an autocatalytic chemical reaction model, that is, the glycolysis model. The fundamental traits of the prescribed numerical method, for instance, the positivity, consistency, stability, etc., are also veriﬁed. The authenticity of the proposed scheme is ensured by comparing it with two extensively used numerical techniques. A numerical example is presented to observe the graphical behavior of the continuous system by constructing the numerical algorithm. The comparison depicts that the projected numerical design is more productive as compared to the other two schemes, as it holds all the important properties of the continuous model.


Introduction
Glycolysis reaction behavior arises in biological sciences, chemistry, biochemical sciences, open chemical reactors, etc. This type of reaction reflects the pattern behavior of the solutions. So, the solution of the underlying reaction model exhibits oscillatory behavior when the steady state of the model is unstable. Hess and Boiteux [1] investigated the metabolic oscillations in various cellular systems. Initially, glycolytic traveling waves were observed in yeast cells. Afterward, the oscillatory wave solutions were observed in many phenomena, such as an open spatial reactor [2] and the glycolysis reaction for the uptake of glucose in cones [3]. The dynamical systems are generally associated with the real-world problems, represented by the mathematical systems involving the partial or ordinary differential equations. Mostly, these differential equations contain the space and time derivatives of the unknown quantities which represent the dynamics of the physical model. In the system of differential equations, the space derivatives of a different order have a strong impact on the behavior of the dynamical system. These types of mathematical models with the time derivatives have vast applications in the various fields of the real world, for instance, chemistry, physics, biological sciences and engineering [4][5][6][7][8][9][10][11][12][13]. Physical chemistry also has a verity of applications in mathematics. There are several chemical reaction models in the form of differential equations that represent the different properties of the chemical substances. Glycolysis is a fundamental chemical reaction which appears in the cytosol of cells' cytoplasm. It is actually a metabolic pathway for the energy of cells in an organism [14,15]. This pathway concerns the oxidative disruption of one glucose molecule into two pyruvate with the addition of energy in the form of adenosine triphosphate (ATP) and nicotinamide adenine dinucleotide (NAD) [16]. Glycolysis plays a significant role in the body cells as glucose is a major source of propellant for the tissues in the body. For instance, the only mean of energy for the brain is the glucose. The body must carry a sufficient quantity of glucose in the blood to guarantee the brain functioning naturally. So, glycolysis appears in almost all living organisms.
To study the advection-diffusion systems, the applications of the derivatives, ultimately the differential equations in the mathematical models of the underlying systems, have a significant role. The governing equations involved in the phenomenon provide an easy understanding for investigating the model. For the numerical solution of the glycolysis model, it is of great significance to verify the important properties such as the existence, uniqueness and convergence toward the true steady states. To explore the numerical solution of the set of equations representing the system, the awareness of the existence and uniqueness of the solution is an important thing. Moreover, the study of the convergence behavior of the solution plays a vital role in the numerical analysis [17]. Before computing an approximate solution, it is more beneficial to ensure the existence of the solution. When the solution is guaranteed, then it is a further challenging task to find how many solutions are possible. Then, a feasible solution is traced, which is in accordance with the constraints of the system. The solutions of differential equations generally lie in the function spaces called Banach spaces. Generally, the desired values of the state variables involved in the system of differential equations cannot have a global bound, so it is obviously attractive to consider a closed ball in the space of continuous functions. The advantage of considering such subsets is that explicit estimates can be obtained. In the first stage, the problems are reduced to the operator having a fixed-point property for which the solutions are guaranteed in the closed balls of these spaces of the functions that lead to the limitations on the radius of the closed balls for an available interval or boundary conditions. Likewise, the same radius chosen in the same problem can also lead to the limitations for admissible auxiliary data. The closed ball for which the limitation or restriction is placed is called the optimal ball. Tutschke, in 2005, gave the idea about explicit estimates for the general equations of operators [18]. Numerical schemes are immensely applied to approximate the solutions of the system of differential equations. So, it is necessarily required to establish such numerical schemes that carry all the fundamental traits of the model, namely the positivity, boundedness and convergence toward the fixed point. In recent years, a number of researchers studied the glycolysis models [19,20] which are nonlinear in nature. Due to the nonlinearity, to find the analytical solutions for these models is a strenuous job, so most of the researchers found the numerical solutions of these types of mathematical models [21,22].
In 1968, Selkov [23] constructed a simple model of glycolysis that represents the qualitative behavior of most experimental data for a single-frequency oscillation in glycolysis. Selkov's model is described as Mickens, in [24], analyzed the above Selkov model numerically and obtained the approximate solution by using a finite difference scheme in a nonstandard way (NSFD), which possesses the property of positivity preserving naturally. In the abovementioned model, u 1 and u 2 represent the concentrations of the chemical substances, so they must be positive. The systems of first-order partial differential equations have been considered within the framework of the Clifford analysis, where the components of the Clifford algebravalued functions are the solutions of the diffusion equations [25]. Verveyko and Verisikin demonstrated the computational analysis of the glycolysis reaction model in an open spatial reactor [26]. In [27], Nauman et al. considered a one-dimensional coupled autocatalytic glycolysis model and solved it numerically by an efficient finite difference explicit numerical scheme, which preserves the positivity of the unknown quantities used in the model. Now, we convert the model, mentioned above, into the advection-reaction-diffusion model as follows: With the inclusion of the advection and diffusion terms in the mathematical model concerning population dynamics and an autocatalytic chemical reaction system, it represents more realistic behavior as advection describes the directed movements and the random movements are represented by the diffusion process. So, the extended models thoroughly explain the relevant physical phenomena [28][29][30][31][32][33][34][35]. In this paper, we dealt with a one-dimensional glycolysis model with advection-diffusion.

∂ f ∂t
∂g ∂t with initial and boundary conditions where f and g are the unknown quantities representing the concentration of chemical substances, and ν 1 and ν 2 represent the diffusion coefficients. In the above proposed model, The remainder of this article consists of the following sections. In Section 2, the general form of the advection-reaction-diffusion system is considered for the analysis. The existence of the optimal solution of this system in an optimal ball of a Banach space with explicit estimates is discussed. A result about the uniqueness of the solution of the proposed system is also established. Section 3 is devoted to the numerical study of the autocatalytic glycolysis model. A reliable and efficient numerical scheme, which is known as a nonstandard finite difference technique, is designed which preserves the traits of the physical system. The consistency and stability of the proposed scheme are also discussed in this section. In order to show the efficiency of the proposed scheme, its comparison with the other well-known numerical schemes is necessary. So, to compare our results, two other existing methods are also applied on the same problem, and the results are compared in this section. A numerical example with complete numerical simulations is illustrated in Section 4. In Section 5, some concluding remarks are given.

Existence of Optimal Solutions
In this section, we find some optimality conditions in Banach spaces for some quantities used in the existence of optimal solutions to glycolysis model (1) and (2). This system can be expressed in a more general form as In the above form, H 1 and H 2 showing in the right sides of Equations (3) and (4) may be linear or nonlinear functions in which nonlinearity of H 1 and H 2 may depend on f and g as well as on their first-and second-order apace derivatives [36]. The equivalent form for the solution of the above system can be written as follows Combining these two equations, we obtain In the light of the theory of operators, the integral equation showing above to evaluate the solution to Equations (3) and (4) can be written in operator form, which becomes a fixed-point problem as follows For any real number τ ≥ 0, let X = C [a, b] × [0, τ] be a Banach space in which we will optimize the proposed problem. Optimization can be performed by choosing a closed and convex subset, called closed ball, in the solution space X. To choose this closed and convex subset in X, two possibilities arise as follows: 1.
Consider a closed ball in the Banach space X whose center is at Θ ∈ X which is usually known as the zero of the space X.
The fixed point results from the literature of a nonlinear analysis can be helpful to ensure the existence of the solution. In this regard, Schauder fixed-point theorem is applied to guarantee the existence result for the system (1) and (2). First, it is important to consider that the fixed-point operator J maps continuously the ball S r (Θ) in to itself.

Now, consider a mapping
Suppose that To fulfill the first condition of Schauder's fixed-point theorem, J is needed to be a self-map. So, If we choose the initial condition f 0 and τ, length of the interval for continuity, we can optimize r, that is, the optimal ball. Then, If we fix r and | f 0 |, then To obtain the highest value of τ, we have to find the maximum value of . Let r * be the value of optimized radius r, then it satisfies the equation (by using (7)) Again, if τ is chosen first, the value of f 0 can be optimized from (6), that is, For the greatest value of f 0 , we have the maximum value of r − τk(r). The optimized value of r satisfies the equation 2.
If we consider an optimal ball with center at the initial condition f 0 .
By Schauder's fixed-point theorem, existence of the solutions to Equations (1) and (2) require the self-mapping first for J : For the largest possible value of τ, maximization technique can be applied from applied calculus on r k 1 (r) in which k 1 (r) is considered as a differentiable function. Suppose that r * be the value of optimal radius, then it should be the root of the equation If If we suppose k 1 (r) is positive for all values of r, then r k 1 (r) has a maximum value at every value r * which satisfies Equation (11). Because between two maximum values there exists at least one minimum value, the assumption of k 1 (r) ≥ 0 implies that Equation (11) has at most one solution.
With the help of the above discussion, we can develop a result.

Theorem 1.
Suppose that for a real-valued continuously differentiable function k 1 (r), the optimal value of radius r * is the root of Equation (11). Moreover, k 1 (r) ∈ C 2 and positive. Then, at most, one value of optimal radius of a closed ball exists in which the solutions of Equations (1) and (2) exist.
Numerical solutions have a key role in the applied study of mathematics. Along with the approximate solutions, numerical simulations help us to understand the behavior of the solutions as well as the numerical scheme [37][38][39]. In the following section, we use a suitable and reliable Mickem's technique, named the nonstandard finite difference (NSFD) scheme, to obtain the solution of the prescribed model, and to show the strong reasonings to use this scheme, we will find the solution to the same problem by some well-known numerical schemes for the comparison.

Numerical Solutions of Glycolysis Model with Analysis
In the literature of numerical analysis, the finite difference schemes have a pivotal role to find the numerical solution of differential equations that may be linear or nonlinear. The use of such kind of techniques makes the numerical calculations very easy as compared to the other existing methods because this technique converts the continuous model into a discrete formulation by taking some finite number of values of the function defined on the corresponding finite number of points in the given domain which is easy to handle [40][41][42][43][44]. The Taylor's series expansion is the best tool to obtain the numerical approximations.
In the remaining part of this paper, assume that m, n are two positive integers and τ ∈ R, a positive real number. For finding the values of the state variables involved in Equations (1) and (2), we proceed as follows: the intervals [a, b] and [0, τ] are discretized as given below a = a+b 2 ). Strogatz discussed the stability of the fixed point of the glycolysis model in [45]. According to Strogatz, the fixed point of glycolysis is stable if a+b 2 < 0. As for some nonlinear systems, it is difficult to investigate the analytical solutions to Equations (1) and (2). So, we adopt some numerical techniques to approximate them. A numerical scheme is called reliable if it sustains the fundamental traits of the system. In particular, it should preserve the structure of the continuous system. Mickens established the criteria for constructing the structure-preserving nonstandard numerical scheme [46]. Here, we apply the scheme which was established by Charpentier [47]. This scheme retains the positivity property for the underlying system. The efficiency of the proposed scheme is elucidated by comparing it with other existing schemes in the literature.
There are many implicit and explicit schemes in the literature proposed by various researchers. In the explicit scheme, the dependent quantities are expressed in terms of some already known quantities for the previous step of time, while for the implicit schemes, both dependent and independent variables (the value of functions at a future time step) are used in one equation, so they can be represented in the matrix notation. We can discretize our proposed model (1) and (2) by using the formulas Substituting these formulas in Equations (1) and (2) and applying the rules defined by Mickens [46], we obtain After some simplifications, we obtain

Order of Accuracy
For the accuracy of the numerical solution, it is very important that is achieved by the numerical design. The accuracy of the proposed scheme is obtained by using Taylor's expansions. From Equation (12), let Similarly, after finding the values of f (x − λ, t + µ), f (x + λ, t + µ) by using Taylor's series, we obtain By applying Taylor's expansion series formulas on (17) and after some simple calculations, we obtain g → ∂g ∂t + ∂g ∂x − ν ∂ 2 g ∂x 2 + ag + f 2 g − b as µ → 0 and λ → 0. Hence, the proposed scheme is consistent with order 1.

Stability of Proposed Scheme
From Equation (14), After linearization, we have which is less than 1.
Again, from Equation (15), after linearization, we obtain Similarly, we obtain which is also less than 1. Hence, the proposed scheme is von Nuemann stable.

Positivity of Proposed Scheme
M-matrix theory [48] is quite useful for proving the positivity of a large number of mathematical models regarding population dynamics, economy, engineering, autocatalytic chemical reactions, etc.
M-matrix over R is a square matrix if all the entries in off-diagonal position are non-positive.

Remark 1.
It is evident that M-matrix obeys the diagonal dominance property and all the entries of the inverse of M-matrix are always positive.
In the presence of all the abovementioned properties, both the matrices K and L are Mmatrices. This implies that K and L are non-singular matrices. So, the Equations (18) and (19) can be represented by The corresponding difference equations are determined by Equations (1) and (2), and after substituting the above approximations, we have

Numerical Experiment
In this section, we discuss a numerical experiment. Consider the system (1) and (2) Figure 1 show the negative solution of g, which is worthless, while this scheme is mathematically stable and unconditionally consistent. Moreover, the graphs in Figure 2 point out that this technique is divergent for various step sizes. So, we can conclude that the Crank Nicolson method is not applicable to the proposed model.
The simulations in Figures 3 and 4 for the upwind implicit method are similar to the Crank Nicolson method. The upwind implicit scheme exhibits the non-physical behavior in these graphs. Therefore, this method also fails to sustain the structure of the glycolysis model.  Graphical solution of f , g (concentration of substances) using upwind implicit method at a = 3.5, b = 0.25, χ 1 = 2, χ 2 = χ * 2 = 0.04 ν 1 = ν 2 = 10 −3 .
The graphical approach of the suggested NSFD implicit technique is performed in Figures 5 and 6. All the important attributes of the state variables involved in the system (1) and (2) are contained by the proposed technique. This technique achieves the convergent solutions for all the step sizes and also holds the positive solution.
In Figures 1-6, we select the values of a and b in such a way that the condition   Graphical solution of f , g (concentration of substances) using designed NSFD implicit method at a = 3.5, b = 0.25, It is evident from the above figures that the underlying system (1) and (2) is unstable a+b 2 is less than zero. Figure 7 validates that the proposed technique demonstrates the same behavior as possessed by the model (1) and (2). In the graphs of Figure 7, we take the values of a and b in such a way that the expression a+b 2 is less than zero. Clearly, the graphical solution exhibits unstable behavior.  Graphical solution of f , g (concentration of substances) using NSFD implicit method at a = 3.5, b = 0.25, χ 1 = 6, χ 2 = χ * 2 = 0.12 ν 1 = ν 2 = 10 −3 .  Graphical solution of f , g (concentration of substances) using designed NSFD implicit method at a = 0.008, b = 0.6, χ 1 = 2, χ 2 = χ * 2 = 0.04 ν 1 = ν 2 = 10 −3 .

Conclusions
In the current article, the existence of the solution for a chemical reaction system, known as a glycolysis model, is described by using a renowned result of the fixed-point theory. Firstly, we considered a function space (solution space) of the continuous functions and optimized the solution by a closed ball which is a subset of the function space. In the framework of optimization, we constructed a result with suitable conditions. The uniqueness of the solution in this regard was also discussed. After acquiring the existence and uniqueness of the values of the state variables, we constructed a scheme for the numerical approximation of the solutions of the system which holds the unconditional positivity of the required solutions. We also proved that our proposed numerical scheme is consistent with the model (1) and (2). The comparison of the current scheme with the existing schemes verified that our proposed scheme is more efficient. The upwind and Crank Nicolson schemes were selected for the comparison. Moreover, these schemes did not preserve the structure of the continuous model. Because the elements of the function space are not globally bounded, it is better to choose the optimal balls for handling the issue. These balls provide the explicit estimates for the approximate solutions to the system. In the future, numerical solutions can be optimized in this manner. Our proposed scheme may be used as a significant tool for solving the different nonlinear systems.