On Nonlinear Reaction-Diffusion Model with Time Delay on Hexagonal Lattice

: In the work, a nonlinear reaction-diffusion model in a class of delayed differential equations on the hexagonal lattice is considered. The system includes a spatial operator of diffusion between hexagonal pixels. The main results deal with the qualitative investigation of the model. The conditions of global asymptotic stability, which are based on the Lyapunov function construction, are obtained. An estimate of the upper bound of time delay, which enables stability, is presented. The numerical study is executed with the help of the bifurcation diagram, phase trajectories, and hexagonal tile portraits. It shows the changes in qualitative behavior with respect to the growth of time delay; namely, starting from the stable focus at small delay values, then through Hopf bifurcation to limit cycles, and ﬁnally, through period doublings to deterministic chaos.


Introduction
Currently, reaction-diffusion models are used in the design and study of many detecting, measuring, and sensing devices. Immunosensors, which are studied here as an example, represent a kind of them. Such spatial-temporal models are described by the systems of partial or lattice differential equations.
The reaction-diffusion models are traditionally studied from the viewpoint of their qualitative analysis. Even in the case of a small number of spatial elements, they show complex behavior. In [1], it was shown that the model describing the chemical reaction of two morphogens (reactants), one of them diffusing within two compartments, resulted in "bi-chaotic" behavior. The origin of such chaotic phenomena (They called it "spiral turbulence" in [2]) were also explained with the help of the statistics of topological defects [2].
When considering continuously-distributed reaction-diffusion models described by nonlinear partial differential equations, the Feigenbaum-Sharkovskii-Magnitskii bifurcation theory can be applied, which results in a subharmonic cascade of bifurcations of stable limit cycles [3].
The lattice differential equations describe systems with a discrete spatial structure, which is more consistent with pixel devices. These equations were also called earlier by a series of authors spatially-discrete differential equations [4].
Due to [5], a typical lattice differential equation takes the form:

Lattice Reaction-Diffusion Model on a Hexagonal Biopixel Array
The model, which is further presented, describes the antigen-antibody interaction, which is the base of the immune reaction used for the construction of immunopixels. Such a model can be applied to general predator-prey models of population dynamics. The model was studied in the case of the rectangular lattice in [17]. Therefore, here, we will compare numerically the influence of the hexagonal lattice on the stability properties of the immunosensor. Therefore, we consider hexagonal lattice of immunopixels, which are indexed by the three-element coordinates (i, j, k), where i, j, k = −N, N, i + j + k = 0, N ∈ N. Let V i,j,k (t) be the antigen concentration, F i,j,k (t) be the antibody concentration within biopixel (i, j, k), and i, j, k = −N, N, i + j + k = 0. Given biopixel (i, j, k), we make the following assumptions with respect to the immune reaction.

2.
We have some probability rate γ > 0 for the binding of antigens with antibodies. 3.
The antigen population tends to some carrying capacity with the rate δ v > 0.

5.
As a result of immune response, we have the increase of the density of antibodies with probability rate ηγ. 6.
The antibody population tends to some carrying capacity with rate δ f > 0.

7.
The immune response appears with some constant time delay τ > 0.
Taking into account the prerequisites mentioned above, we obtain a simplified antibody-antigen competition model with delay for a hexagonal array of biopixels, which uses the Marchuk model of the immune response [18][19][20][21] and the spatial operatorŜ, which is constructed similarly to [22] (Supplementary Information, p. 10).
with given initial functions: (3) Figure 1. Diffusion of antigens for the hexagonal lattice model. Antigens from six neighboring pixels interact. n > 0 is the constant of disbalance. Here, '1', '3', '5', '8', '9', '11' have to be replaced We use the following spatial operator of discrete diffusion for a hexagonal array of pixels: Each pixel is affected by the antigens flowing out of six neighboring pixels, two in each of three directions of the hexagonal array. The adjoint pixels are separated by the distance ∆.
Boundary conditions V i,j,k = 0 for the edges of the hexagonal array, i.e., if i ∨ j ∨ k ∈ {−N − 1, N + 1}, are used. Sections 3 and 4 present analytical results with respect to the model (2) in the form of restrictions for the parameters, enabling us persistence and global asymptotic stability. Moreover, Section 5 offers numerical research of the system's qualitative behavior with dependence on the changes of the time of immune response τ (time delay), diffusion rate D∆ −2 , and factor n.

Persistence of the Solutions
We use the following definition, which generalizes [23] for lattice differential equations. Definition 1. System (2) is said to be uniformly persistent if for all i, j, k = −N, N, i + j + k = 0, there exist compact regions D i,j,k ⊂ int R 2 such that every solution (V i,j,k (t), F i,j,k (t)), i, j, k = −N, N, i + j + k = 0 of (2) with the initial conditions (3) eventually enters and remains in the region D i,j,k . then: for some large values of t. Here: Proof. Firstly, we can prove that there exists some large instant of time Let us assume the contrary, i.e., there is i , j , k ∈ −N, N, i + j + k = 0 such that S{V i ,j ,k (t)} > 0 at t > T 1 , which is a contradiction with the balance principle.
Since the solutions of the system (2) and (3) are positive, then: Further, we can apply the basic steps of the proof of Lemma 3.1 [24], which is proven in the non-lattice case (i.e., without the spatial operator).

Remark 1.
The conditions of the uniform persistence of System (2) in the non-lattice case were obtained in [25]. They resulted in Inequality (5) provided that: holds.

Stability Problem in Immunosensors
Stability from a biosensor point of view is very important. With respect to biosensors, two notions of stability are introduced. Namely, these are self-stability and operational stability.
Self-stability characterizes the sustainability of the biosensor for a long period of time. To check the self-stability of the biorecognition element, measurements of the response of the biosensor to the same amount of measurands are fulfilled.
Operational stability shows the identity of the sensor response to the same substrate concentration when conducting a large number of consecutive measurements [26] and is closely related to the metrological characteristic of convergence (repeatability). The operational stability of the sensors is characterized by a relative standard deviation with repeated measurements of the standard sample and the residual activity of the biosensor after several consecutive measurements over a certain period of time.
The results, which are presented in the work, could be applied to study both types of biosensor stabilities. Namely, with help of the changing initial conditions, we can simulate different circumstances to investigate its qualitative behavior.

Steady States
The steady states can be calculated as a result of solving dV i,j, Prey-free steady state: In case of the absence of antigens (preys), The last one is of biological meaning and is not reachable for positive initial values (3). In order to get an endemic state (2), the following system of algebraic equations has to be solved: The solving of (10) with respect to V * i,j,k , F * i,j,k results in the system of algebraic lattice equations for V * i,j and taking into account the dependence Analyzing the solutions, two types of endemic steady states are differed. Pixel-independent endemic state: In the case of the existence of solutions of (10), which are identical for all pixels, i.e., provided that ηγβ > δ v µ f . Pixel-dependent endemic state: Because of the diffusion between pixels D/∆ −2 , the endemic steady states are not equal to (11), and in Section 5, it is evidenced numerically.
For the diffusion-free system, i.e., D = 0, pixel-independent endemic states are observed for pixels on the hexagonal boundary.
If D > 0, we get pixel-dependent endemic steady states for the internal pixels, which converge to the pixel-independent ones (11). Numerical study shows that the bigger the number N, the clearer is the appearance of pixel-dependent endemic states.

Global Asymptotic Stability
Here, the global asymptotic stability of the positive endemic steady state E * i,j,k = V * i,j,k , F * i,j,k , i, j, k = −N, N, i + j + k = 0 is investigated. Theorem 1. Assume that: and: Then, the positive endemic steady state E i,j = V i,j , F i,j , i, j = 1, N of system (2) is globally asymptotically stable.

Proof.
We use the change of the variables: This implies: Substituting (13) into (2), we get: This yields: Since: then: Rewrite the first equation as: For arbitrary pixel (i, j, k), i, j, k = −N, N, i + j + k = 0 consider the Lyapunov functionals: Then, we have for the right-upper derivative along (18): By Lemma 1, we know that there exists a T 0 ≥ 0 such that V(t) ≤ M v for t ≥ T 0 , and hence, for t ≥ T = T 0 + τ, Next, let: Then, we have at t ≥ T: Rewrite the second equation (17) as: For arbitrary pixel (i, j, k), i, j, k = −N, N, i + j + k = 0, consider the Lyapunov functionals: Then, we have for right-upper derivative along (20): Considering for arbitrary pixel (i, j, k), i, j, k = −N, N, i + j + k = 0 the Lyapunov functional G i,j,k (t) = G i,j,k,1 (t) + G i,j,k,2 (t), we have for the right-upper derivative along (20): Let us consider: Since V i,j,k is the steady state, we see that for any arbitrary small > 0, there is some instant of time T 1 ( ) such that:S Further, we consider t ≥ max {T 0 , T 1 ( )} and omit in (23) without lost of generality. We have: Then, calculating the derivative of G(t) along the solutions of (20), we have: |e y i,j,k (t) − 1| ≤ 0.

Remark 2.
The conditions of the global asymptotic stability of the endemic steady state E i,j,k = V i,j,k , F i,j,k , i, j, k = −N, N, i + j + k = 0 of system (2) are dependent on the diffusion rate D, the distance between pixels ∆, and the disbalance constant n.

Numerical Study
For the numerical simulation, we consider Model (2) of the hexagonal pixel array at N = 4, β = 2 min −1 , γ = 2 mL min·µg , µ f = 1 min −1 , η = 0.8/γ, δ v = 0.5 mL min·µg , δ f = 0.5 mL min·µg , D = 0.2 nm 2 min , ∆ = 0.3 nm. Numerical modeling was implemented at different values of n ∈ (0, 1]. For this purpose, we used the RStudio environment. Using the local bifurcation plot, the dynamics of the system (2) was analyzed for different values of n ∈ (0, 1] (see Figure 2). We concluded that oscillatory and then chaotic behavior starts for smaller values of τ at smaller values of n. Further, increasing the values of n, we can observe asymptotically stable steady solutions for a wider range of τ. Numerical integration of the system have shown the influence of time delay τ. Namely, as it agrees with our analytical results of Section 4, we observe the stable focuses at pixel-dependent endemic states for small delays τ ∈ [0, 0.18) (see Figure 3). At τ ≈ 0.18 min, the stable focus is transformed into a stable limit cycle of tiny radius, which corresponds to Hopf bifurcation. A deeper study of this phenomenon requires obtaining the condition of the appearance of the pair of purely imaginary roots of the characteristic quasipolynomial of the linearized system. The limit cycles of the ellipsoidal form are observed till τ ≈ 0.285 min (see Figure 4 for τ = 0.25). Notice that when increasing τ, near τ = 0.285, we get period doubling (see Figure 5) (It can be approximately seen from local bifurcation plot also.).  The qualitative behavior of the immunosensor model can be analyzed with the help of hexagonal tiling plots as well. For this purpose, we can use both plots for antigens ( Figure 6), antibodies (Figure 7), and the probabilities of binding antigens by antibodies (Figure 8).

Conclusions
In this work, a reaction-diffusion model of a hexagonal immunosensor was considered. Mathematically, it was described by the system of lattice delayed differential equations on the hexagonal grid. The system included the spatial operator describing the diffusion of antigens between pixels.
The main results dealt with the qualitative investigation of the model. The conditions of the global asymptotic stability were obtained. These used the construction of the Lyapunov functional. These resulted in an inequality, including the system parameters and delay. Therefore, the estimation of the delay enabling the global asymptotic stability can be obtained.
Numerical analysis of the model was performed with the help of the bifurcation diagram, phase trajectories, and hexagonal tile portraits. The changes in the qualitative behavior with respect to the growth of the time delay were shown; namely, starting from the stable focus at small delay values, then through Hopf the bifurcation to limit cycles, and finally, through period doublings to deterministic chaos. This agreed with the results on space-time chaos for reaction-diffusion systems, which were previously obtained in [1][2][3].
As compared with the rectangular lattice model studied in [17], for the hexagonal model, we observed Hopf bifurcation at smaller values of τ. That is, the hexagonal lattice accelerated changes in qualitative behavior.
Note that the model can be applied for an arbitrary amount of pixels determined by natural N ≥ 1. However, it can be numerically seen that the qualitative behavior of the entire immunosensor was determined by a seven-pixel array.